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Abstract. In these lecture notes, we present a pedagogical review of a number of related numerically exact approaches to 
quantum many-body problems. In particular, we focus on methods based on the exact diagonalization of the Hamiltonian 
matrix and on methods extending exact diagonalization using renormalization group ideas, i.e., Wilson's Numerical Renor- 
malization Group (NRG) and White's Density Matrix Renormalization Group (DMRG). These methods are standard tools 
for the investigation of a variety of interacting quantum systems, especially low-dimensional quantum lattice models. We also 
survey extensions to the methods to calculate properties such as dynamical quantities and behavior at finite temperature, and 
discuss generalizations of the DMRG method to a wider variety of systems, such as classical models and quantum chemical 
problems. Finally, we briefly review some recent developments for obtaining a more general formulation of the DMRG in 
the context of matrix product states as well as recent progress in calculating the time evolution of quantum systems using the 
DMRG and the relationship of the foundations of the method with quantum information theory. 



1. INTRODUCTION 
1.1. Strongly Correlated Quantum Systems 

Many problems in condensed matter physics can be described effectively within a single-particle picture |1]. 
However, in systems where interaction leads to strong correlation effects between the system's component parts, 
mapping to an effective single particle model can lead to an erroneous description of the system's behavior. In such 
cases, it is necessary to treat the full many body problem. As we will see later, it is generally impossible to treat the full 
ab initio problem, even numerically. However, in the investigation of condensed matter systems one is often interested 
in the low-energy properties of a system. In this article, we will discuss numerically exact methods to treat effective 
models that describe these low-energy properties. 

One possibility to numerically investigate a physical system is to discretize the underlying differential equation or to 
map the model onto or formulate it directly on a lattice. It is then necessary to describe model properties on the lattice 
sites (e.g., is it a spin system? Are the particles fermions or bosons?) and to specify the topology of the lattice, its 
dimension and the boundary conditions one needs to apply (e.g., a one-dimensional chain or ring, a two-dimensional 

square lattice, a honeycomb lattice, etc.). The system is then composed of N quantum mechanical subsystems, where 

(£) 

each subsystem is located on a site j and is described by a (usually finite) number of basis states |ce- ), £ = 1, . . . 
which depend on the model and may vary from site to site. This basis will be referred to as the site basis in the 
following and is usually "put in by hand" in order to model a physical situation or system. 

One often is interested in general properties of the physical models that go beyond their behavior on finite lattices. 
In these cases, it is necessary to investigate the behavior in particular limits, e.g., by investigating the model's behavior 
when S( — > °° (the continuum or thermodynamic limit) or N — > °° (the thermodynamic limit). Insight can also be 
gained by considering the limit of a continuous quantum field by taking I — > x/a, where a is the lattice constant. In 
the following, we will restrict ourselves to the description of quantum systems on a finite lattice with N sites, as, e.g., 

realized in solid state systems or in recent experiments on optical lattices Q|. Given the site-basis {|cty )}, a possible 
configuration 1 0Cj , oif 1 otffl ) of the total system is given by the direct product of the component basis, 
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In the following, we will omit the state-index I. The set of all of these states constitutes a complete basis and hence 
can be used to represent an arbitrary state 



l*> = E Y(ai,a2,---a N )\a l ,a 2 ...a N ) (2) 
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where the total number of possible combinations and therefore the dimension of the resulting Hamiltonian matrix 
or state vector is Yl 1 j=i s j- This exponentially increasing dimension of the basis of the system is the key problem to 
overcome when treating quantum many-body systems. In principle, a quantum computer would be needed to provide 
the full solution of an arbitrary quantum many -body problem L3j|4j|5|]. However, a variety of numerical approaches 
exist to investigate such systems exactly on ordinary computers. In this article, we restrict ourselves to a set of related 
techniques; specifically, we will discus exact diagonalization and numerical renormalization group (NRG) methods, 
including both Wilson's original NRG |6] and the density matrix renormalization group (DMRG) |7|,|8|]. 

The authors are participants in the ALPS collaboration, whose aim it is to develop and provide free, efficient, and 
flexible software for the simulation of quantum many body systems f|f]. In the ALPS project, implementations of a 
number of different numerical methods for many-body systems are available - including exact diagonalization and 
DMRG. A variety of models on different lattices can be treated; the programs themselves can be used as a "black 
box" to carry out calculations or can be a helpful starting point for developing one's own implementation. The reader 
is encouraged to visit the URL and to try out the programs. 



The behavior of a quantum mechanical system is, in general, governed by the time-dependent 



or time-independent Schrodinger equation, 



ih?- t \V(t))=H\V(t)) (3) 



H\V>) ^E^} . (4) 



One approach to treating quantum mechanical systems numerically is to solve (at least partially) the eigenvalue 
problem formed by the system's Hamiltonian. Note that Hamiltonians without interaction, i.e., consisting of a sum 
of single-particle terms, can be solved by treating the single-particle problem - there is no need to work in the full 
many body basis. The many-body wave function can then be constructed out of the single-particle solutions. However, 
for some systems, e.g., systems with strong interactions, it is not possible to reduce the problem to a single particle 
problem without introducing substantial and, in many cases, uncontrolled approximations. The approaches discussed 
in this article are suitable for such systems and have been applied successfully (mostly for one-dimensional systems) 
to investigate the low-energy behavior. 

A typical Hamiltonian for a quantum lattice system can contain a variety of terms which can connect an arbitrary 
number of subsystems with each other. A general form for quantum lattice Hamiltonians is given by 

H = £HW+£fl2 ) + -+ I <, + ...- (5) 
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The local or site Hamiltonian Hj 1 ^ describes general, local model properties of the system and may contain additional 

(2) 

terms such as a local chemical potential or coupling to a local external field. The two-site terms H [m connect pairs 
of subsystems. Such terms can connect two sites or, in the case of impurity systems or multi-band systems may also 

(2) 

contain hybridization or couplings between orbitals on the same site. The coupling associated with Hi' is often short- 
ranged (between nearest neighbors or next-nearest neighbors only) in models of interest. Hamiltonians containing only 
terms of the type Hi separate into single-particle problems. In some cases, a change of the single-particle basis can 
transform terms of the type into terms of the type flj , e.g., the diagonalization of a hopping term by Fourier 

(4) 

transformation. However, two-site terms in general lead to an interacting many-body system. Four-site terms Hf mnp 
generally cannot be reduced to a sum over single-site terms; as demonstrated in the following, they always couple at 
least two particles and therefore preclude separation into a single particle basis. The interaction terms of the strongly 
correlated quantum systems treated here will contain terms of the type and sometimes H^ np . 



AAAA 



FIGURE 1. Left-hand side: sketch of the pseudo-potential for a chain of atoms and the corresponding localized Wannier orbital 
centered on the positions of the atoms, which are assumed to be fixed (Born-Oppenheimer approximation). Right-hand side: sketch 
of the effective lattice model. Since the overlap between the Wannier-orbitals falls off very rapidly with distance, it is often sufficient 
to include only the nearest-neighbor hopping t, which leads to the simple tight-binding-model 1 1], 



1.1.1. Example Hamiltonians 

Roughly speaking, one can classify the models of interest as fermionic or bosonic lattice models, as pure spin- 
models, and as impurity models. To model a regular solid, one usually works within the Born-Oppenheimer approx- 
imation, i.e., one assumes the atomic nuclei to be moving slowly enough so that they can be considered to be fixed 
on the relevant time scale for the electronic problem. The simplest fermionic lattice model is the tight-binding model 
(jl, which treats electrons as being localized to Wannier-orbitals centered on a regular array of sites. Such a model 
would be appropriate for unfilled d or / orbitals in transition metals, for example. In the simplest case, one restricts 
the description to only a single band, although multi-band models can also be considered. Due to the finite overlap 
between nearby orbitals, the particles can "hop", from one site to another, i.e., tunnel between the corresponding Wan- 
nier orbitals, with an amplitude t . The mapping of unfilled strongly localized orbitals to the tight-binding model is 
sketched in Fig.^ 

Since the overlap between strongly localized orbitals typically decreases exponentially with separation, in most 
cases one considers hopping only between nearest-neighbor sites, or at most hopping to next-nearest neighbor sites. 
The Hamiltonian for the simple tight-binding chain with nearest-neighbour hopping in second quantization is given by 

^ tb = -L^(4 + i c ; + c >;+i) • (6) 

(hi) 

Here and throughout this article, (i, j) will denote pairs of nearest-neighbor sites, and we will choose units so that h=\. 
The quantity fy is the tunneling matrix element between the sites i and j and corresponds to the kinetic energy of the 
particle "hopping" between these two sites. The operators cj (c .) create (annihilate) fermions in the subsystem labeled 

by j, and obey the usual anticommutation relations, {cj ,cj} — 8jj and {cj,cj} = {c,-,c ; } = 0. Without additional 
interaction terms, this model can be reduced to a single -particle problem, which will be treated in Sec. 13.41 

When spin is included, Pauli's exclusion principle allows up to two fermions with opposite spin to occupy a Wannier 

orbital. The possible states on a site j are then {|aj^)} = {|0/), | |y), | J.y), KT4);)}- This is an appropriate site basis 
for modeling systems with Coulomb interaction, e.g., between electrons, 

9 - e2 

Hjm ~\rj-r m \- (?) 

One prominent example is the Hubbard model, 

H=-t £ {c) a c ma +cj na c jc ^ +UY j n j ^n Jil , (8) 

U,m).o j 

with rj labeling the spin and nj. a = Cj a cj. a the local particle number operator. Since screening effects are strong in 
strongly localized d and / orbitals, it is usual to reduce the Coulomb interaction to the on-site term (the Hubbard 
term), or sometimes the nearest-neighbor term V £«,«,+,•, in which case the model is called the extended Hubbard 
model. For such models, the dimension of the Hilbert space grows as A N , where N is the number of lattice sites. Since 
the interaction terms couple two particles, it is not possible to reformulate the problem in a single particle basis - an 
exact treatment must treat the full many -body basis | Of i , 0S2 , , <Xn)- The eigenstates of such systems are many -body 



FIGURE 2. The second-order virtual exchange process in the Hubbard model which leads to the antiferromagnetic Heisenberg 
coupling, see Ref. 1 10]. 



states, and the excited states cannot, in general, be mapped to single-particle excitations. Within this formalism, it is 
also possible to model disorder, using, for example, the local Anderson Hamiltonian 

II? ^A,» /0 . (9) 

a 

where A ; - varies randomly from site to site according to a given probability distribution. 

Another important class of models are spin models, where S ( - describes localized quantum mechanical spins (S = 
1/2, 1,3/2, .. .). The possible states on a site j are then (omitting the site index) {|oc^)} = {| — S), \ — S+ 1), . . ., \S)}, 
i.e., such models possess (25 + 1 ) N degrees of freedom. A prominent example is the Heisenberg Hamiltonian, 

H Heis = J £ Sj ■ S m = P £ S)S; n + £ (Sp- + SJS+) . (10) 

0» {;•'«> 0» 

The 5=1/2 Heisenberg model can be derived as the strong coupling limit of the Hubbard model (U /t — > °°) at half 
filling, i.e., average particle density (n) = 1. The corresponding antiferromagnetic exchange coupling in second-order 
perturbation theory (see Fig.|2j is J = 4t 2 /U | 10]. 

Several variations of this model are known. For example, when J xy = 0, it is termed Heisenberg model with an Ising 
anisotropy, or, more generally, with XY anisotropy, if J z ^ J xy . Additional terms such as H" = D(S z j) 2 ("single-ion 

anisotropy"), or a biquadratic coupling, H b ^ n — J2 (Sj ■ S„,) 2 can be introduced. When doped below half-filling, the 
strong-coupling limit of the Hubbard model is described by the t-J model 

H tJ = &H tb &+J £ (Sj ■ S m - jnj n m ) , (11) 

o» ^ 4 / 

where & projects onto the space of singly occupied sites, see, e.g., Ref. 1 10], leading to a site-basis with three states 
per site. This Hamiltonian on a two-dimensional square lattice has been proposed as an effective model for the Cu02 
planes in the high-r f superconductors. 

A third important class of models are Hamiltonians describing systems containing impurities. A prominent example 
is the Anderson impurity model, which models the hybridization of the localized d (or f) orbital of an impurity located 
at a site j with conduction-band orbitals (in the simplest case) located at that site, while the electrons on the impurity 
interact with each other via a Hubbard term. The localized part of the Hamiltonian is then 

Hf = Z £ d n M + V L Ko<-,.o • + '-'<•<. , (12) 

while the conduction band is modeled as a tight-binding band 0. Here dj a is the annihilation (creation) operator 

for a fermion with spin a on the impurity site, and nj a = dj a dj, a . This problem has been treated with one or two 
impurities, leading to the the single or two impurity Anderson model, (SIAM or TIAM, respectively) or as a system 
with impurities located at every lattice site, leading to the periodic Anderson model (RAM). 

Perhaps the most well-known example of an impurity system is the Kondo model, which describes the localized 
impurity level as a quantum mechanical spin S, leading to the on-site Hamiltonian 

H f = Y s J-( c l« a °,P c j,p) (13) 

where a„ R is a vector whose components are the Pauli matrices. This model can be derived as the limit of the 

J — L 

symmetric SIAM at strong coupling, U > (t,V) IllH . 



FIGURE 3. Examples of two-dimensional lattices. On the left and in the middle the square and the Kagome lattices, respectively, 
are shown. The sketch on the right shows how a lattice can be tilted along the symmetry axis - in this way, it is possible to treat 
larger system sizes by exploiting possible additional translation symmetries. 



7.7.2. Lattices 

The lattice structure of a regular solid is described by a Bravais lattice, which is built up from a unit cell and a set 
of translation vectors which create the regular pattern. Since the memory of a computer is finite, it is often necessary 
to restrict the treatment to finite lattices. In this case, boundary conditions must be specified. For open boundary 
conditions (OBC, sometimes they are also calledfixed BC) the intersite terms of the Hamiltonian are simply set to zero 
at the boundaries. This leads to a lattice which is not translationally invariant. Translational invariance can be restored 
on a finite lattice by applying periodic boundary conditions (PBC, also sometimes called Born-von-Karman BC), in 
which the edges of the lattice are connected together in a ring geometry for one-dimensional systems or a torus for 
two-dimensional systems. It is also possible to introduce a phase at the boundaries, leading to antiperiodic or twisted 
boundary conditions. In Fig.[5]some examples of two-dimensional lattices are shown. 

Lattices possess symmetries which can be exploited in calculations. The translation symmetries, parameterized by 
multiples of the Bravais lattice vector, can be used for boundary conditions that do not break translational invariance, 
i.e., (A)PBC or twisted boundary conditions. Other useful and simple symmetries are symmetries under a discrete 
lattice rotation, e.g., rotation by n/2 for a square lattice (group C4,,), or reflections about a symmetry axis. 

In order to have more sizes available for two-dimensional translationally invariant lattices, it is helpful to work 
in geometries where the clusters are tilted with respect to some symmetry axis, as depicted in Fig. [5] For the two 
dimensional square lattices, the spanning vectors for such tilted lattices are given by F\ = (n,m), F2 = (— m,n) and the 
total number of sites is then N = n 2 + m 2 . The lattice periodicity then appears for translations at appropriate values of 
(n,m). Reflection and discrete rotation symmetries are still present in these lattices, although they become somewhat 
more complicated. For methods which treat the entire Hilbert space, such as exact diagonalization, it can be important 
to use these symmetries to reduce as much as possible the dimension of the basis that must be treated numerically. 

2. EXACT DIAGONALIZATION 

The expression "Exact Diagonalization" is used to describe a number of different approaches which yield numerically 
exact results for a finite lattice system by directly diagonalizing the matrix representation of the system's Hamiltonian 
in an appropriate many-particle basis. The simplest, and the most time- and memory- consuming approach is the 
complete diagonalization of the matrix which enables one to calculate all desired properties. However, as shown in 
the introduction, the dimension of the basis for a strongly interacting quantum system grows exponentially with the 
system size, so that it is impossible to treat systems with more than a few sites. If only properties of low- or high-lying 
eigenstates are required, (in the investigation of condensed matter systems one is often interested in the low-energy 
properties), it is possible to reach substantially larger system sizes using iterative diagonalization procedures, which 
also yield results to almost machine precision in most cases. The iterative diagonalization methods allow for the 
calculation of ground state properties and (with some extra effort) some low-lying excited states are also accessible. In 
addition, it is possible to calculate dynamical properties (e.g., spectral functions, time-evolution) as well as behavior 
at finite temperature. Nearly every system and observable can be calculated in principle, although the convergence 



properties may depend on the system under investigation. In the following, we describe the most useful methods and 
their variants, which should allow for the investigation of most systems of interest. 

We also discuss the use of the symmetries of a system, which can be important in reaching the largest possible 
system sizes, because the dimension of the matrices to be diagonalized is then significantly reduced. An additional 
advantage is that the results obtained are resolved according to the quantum numbers associated with the symmetries. 
Useful introductions and discussion of the use of symmetries for specific example Hamiltonians can be found in 
Refs. E2E0- A more general mathematical description in terms of group theory is presented in Ref. 1 15]. 

With present-day computers, system sizes can be treated which are large enough to provide insight into the physics 
of many systems of interest. It is possible to treat 5=1/2 spin models with up to about N = 40 sites; the maximum 
sizes reached thus far are, N = 40 sites on a square lattice, N = 39 sites on a triangular lattice, and N = 42 sites on a 
star lattice geometry. The t-J model on a checkerboard and on a square lattice with 2 holes have been treated on up to 
N = 32 sites. Hubbard models at half filling on a square lattice with up to N = 20 sites have been diagonalized; the 
same size was also reached for quantum dot structures. Models with phonon degrees of freedom such as the Holstein 
model, are much harder to treat because the phonons are bosons, which, in principle, have an infinite number of 
degrees of freedom. By truncating the number of phonon states, it was possible to treat a chain with N =14 + phonon 
pseudo-sites. For these calculations, it was necessary to store between 1.5 — 30 • 10 9 basis states 



2.1. Representation of Many Body States 

In order to represent the many-body basis on a computer, it is convenient to map the states of the site basis to a 

single bit or to a set of bits. Thus, every basis state | af^ , a£p , . . . 0$ ) can be mapped to a sequence of bits. In C or 
C++, the most efficient way to handle a bit sequence is to work directly on the bit representation of a (long) unsigned 
integer variable using bit-operators. 

There is usually a natural bit representation for a particular choice of the site basis, The terms of the Hamiltonian 
then can be implemented as sequences of bit operations on the chosen bit set. For example, a basis state of a Heisenberg 
spin- 1/2 lattice may be mapped to the bit sequence 

I T1I2I3 • • • U-iU) — » I1O2O3 . . . ltf-ii* • 

In C or C++, a spin flip operation can be implemented compactly and efficiently as a sequence of boolean operations 
that invert a particular bit. 

For Hubbard-like models, one natural choice for the site basis to the mapping is 

\n]nI) — >n}nI 

with N° = {0, 1}, i.e., one needs two bits for the mapping of a single site. Alternately, a mapping \NfSl) — ► NfSj 
may also be used (mapping Sj = {±1/2} to the bits {0, 1}). 

Other models can be implemented similarly; for bosonic systems, however, it is necessary to introduce a cut-off in 
the number of particles per site so that the site basis has a finite number of states. Typically, one only keeps a very 
small number of on-site bosons, which may depend on the desired system sizes and the parameter values. A minor 
complication occurs when the number of bits needed to represent a state exceeds the number of bits in an integer; then 
more than one integer variable must be used for the description of a single basis state. 

In order to minimize the dimension of the Hamiltonian matrix, it is essential to exploit symmetries of the system, as 
described in more detail in Ref. 11511 . Given a symmetry group & with generators g p , if 

[H,g P }=0, (14) 

the Hilbert space can be partitioned into sectors corresponding to irreducible representations of the symmetry group 
so that the Hamiltonian matrix H becomes block diagonal. The solution of the eigenvalue problem for each block then 
yields the portion of the spectrum of the system associated with a particular conserved quantum number. 

It is possible to exploit continuous symmetries, such as the conservation of the particle number or the conservation or 
the z projection of the spin, S z . For these Abelian U(l)-symmetries, the conserved quantum numbers correspond to the 
sum of all bits representing a basis state. Therefore, all possible basis states preserving this symmetry can be obtained 
by calculating all possible permutations of the bits representing a suitable basis state. Another important symmetry is 



TABLE 1. The reduction of the dimension of the 
sector of the Hamiltonian containing the ground state 
for the 5=1/2 Heisenberg model on a tilted square 
lattice with \/40 x \/40 sites (from Ref. fl6ll ). 



full Hilbert space: 


dim= 


= 2 40 = 10 12 


constrain to S : = 0: 


dim= 


= 138 x 10 9 


using spin inversion: 


dim= 


= 69 x 10 9 


utilizing all 40 translations: 


dim= 


= 1.7 x 10 9 


using all 4 rotations: 


dim= 


= 430,909,650 



the conservation of the total spin, which is an SU(2) symmetry. This symmetry is more difficult to implement than U(l) 
symmetries because it is non-Abelian. However, the Z2 symmetry corresponding to spin inversion can be used instead 
and is easily implemented. Space group symmetries include translational invariance, which is an Abelian symmetry, 
or point group symmetries such as reflections or rotations, which are non-Abelian in general. A set of representative 
basis states can be formed from symmetrized linear combinations of the original basis states generated by applying 
the appropriate generators flUl . 

As an example of how symmetries can reduce the dimension of the largest block of the Hamiltonian that must be 
diagonalized, we show the size of the sector containing the ground state in Table^for a S =1/2 Heisenberg model on 
a V40 x \/40 cluster. In this case, the dimension of the sector to be diagonalized is reduced by a factor of more than 
2500. 

After considering the symmetries, a set of basis states {|«)} is stored in an appropriate form such as a set of bit 
strings or linear combinations of bit strings. The multiplication is then carried out as 

ff|V>=J>|y>tf|n) (15) 

n 

where the set of coefficients (n\\j/) is a stored as vector of real (or complex) numbers with dimension equal to that of 
the targeted block of H. When the Hamiltonian is applied to a basis state \n), the result is a linear combination of basis 
states 

H\n) =Y,{n'\H\n)\n') , (16) 

n> 

where there are typically only few nonzero terms for a short-ranged Hamiltonian. When the states \n) can be 
represented as bit strings, it is usually easy to identify the bit string corresponding to the {|n')} and to determine 
the coefficients (n'\H\n). However, the bit strings then have to be mapped back to an index into the vector of basis 
coefficients in order to calculate the coefficients {n\\\f). A simple way to implement the needed bookkeeping is to 
index all the basis states \n) by an integer corresponding to the bit string and to store this index in an additional vector. 
In this way, it is easy to identify which element of the coefficient vector has to be modified. This implementation is 
simple and fast, but the index vector is of the dimension of the full bit string used to encode the states and is thus not 
memory efficient. More memory-efficient implementations are possible; for example, the bit strings could be stored in 
a lookup-table along with the associated indices and the lookup could be performed using a hash table 111 711 - 

2.2. Complete Diagonalization 

The diagonalization of real symmetric or complex Hermitian matrices is a problem often encountered in numerical 
methods in physics. For example, as we will see later, in one step of the DMRG procedure it is necessary to diagonalize 
the density matrix, represented as a dense, real, symmetric matrix. A number of software libraries provide complete 
diagonalization routines which take a matrix as input and return all of the eigenvalues and eigenvectors as output. 
Among the m ost p rominent are the NAG library 1 181], the routines published in Numerical Recipes as well as 

the LAPACK 1 20] library, which in combination with the BLAS 12lLl22ll23ll provides a very efficient non-commercial 
implementation of linear algebra tools. Such routines could, in principle, also be used to diagonalize the Hamiltonian 
matrix of a finite quantum lattice system directly. However, they are, in general, substantially less efficient than the 
iterative methods described in Sec. l2.3l for finding the low-lying states of the sparse matrices found for such systems. 

The approach normally used Ill7ll is first to transform the matrix to tridiagonal form using a sequence of Householder 
transformations and then to diagonalize the resulting tridiagonal matrix T using the QL or QR algorithm, which carries 



out a factorization T = QL with Q an orthogonal and L a lower triangular matrix. The computational cost of this 
combined approach scales as l« 3 if only the eigenvalues are obtained, and w 3n 3 when the eigenvectors are also 
calculated, where n is the dimension of the matrix. 

The limitations of using such complete diagonalization algorithms to diagonalize the Hamiltonian matrix are 
obvious: the entire matrix has to be stored and diagonalized. Since the dimension of the Hamiltonian matrix grows 
exponentially with system size even when all symmetries are taken into account, the largest attainable lattice sizes for 
strongly correlated quantum systems are generally very small. For a Hubbard chain, one may reach less than 10 sites 
on a supercomputer, while with the iterative diagonalization procedures (presented in Sec. 12. 3i around 20 sites can 
be reached when most of the symmetries of the system are taken into account. As we will see later, one-dimensional 
systems containing more than 1000 sites can be treated on a standard desktop PC using the DMRG, a method carrying 
out an iterative diagonalization in a reduced Hilbert space. A complete diagonalization of the Hamiltonian matrix is 
nevertheless useful for testing purposes, and if a significant fraction of all eigenstates is needed on small systems. 



2.3. Iterative Diagonalization: the Lanczos and the Davidson Algorithm 

If only the ground state and the low-lying excited states of a system are required, powerful iterative diagonalization 
procedures exist which can handle matrix representations of the Hamiltonian with a dimension substantially (up to 
four or five orders of magnitude for short-range quantum lattice models) larger than complete diagonalization. These 
methods can be extended to investigate dynamical properties, time evolution, and the finite-temperature behavior of 
the system. In addition, they form a key part of the DMRG algorithm, which carries out an iterative diagonalization in 
an optimal, self-consistently generated reduced basis for a system. 

The basic common idea of the different iterative diagonalization algorithms is to project the matrix to be treated H 
onto a subspace of dimension M (where N is the dimension of the Hilbert space in which the diagonalization 
is carried out) which is cleverly chosen so that the extremal eigenstates within the subspace converge very quickly 
with M to the extremal eigenstates of the system. The main approach used in physics is the Lanczos method, while in 
quantum-chemistry the Davidson or its generalization the Jacobi-Davidson algorithm is more widely used. 

A number of standard textbooks treat numerical methods for matrix computations in detail. In this section, we will 
discuss only the basic Lanczos and Davidson methods for handling hermitian eigenvalue problems for pedagogical 
reasons. For further details and for the applicability of the methods to non-hermitian and to generalized eigenvalue 
problems, we refer the reader to the following sources: Refs. Q |25|,|26] discuss numerical methods for matrix 
computations in general and also include iterative diagonalization algorithms, including the more modern Jacobi- 
Davidson algorithm (at least in the later editions). A nice overview and a compact representation of the algorithms can 
be found in Ref. [27]. The mathematical theory of the Lanczos algorithm has been worked out in Ref. |28|. 

The basic idea behind iterative diagonalization procedures is illustrated by the very simple power method. In this 
approach, the eigenpair with the extremal eigenvalue is obtained by repeatedly applying the Hamiltonian to a random 
initial state |vo), 

\v„)=H n \v )- (17) 

Expanding in the eigenbasis H\i) = Xi\i) yields 

|v„) = Z,{i\v )H n \i) 

i 
i 

It is clear that the state with the eigenvalue with the largest absolute value will have the highest weight after many 
iterations n, provided that |vo) has a finite overlap with this state. The convergence behavior is determined by the 
spacing between the extremal eigenvalue and the next one; the contribution of the state with the eigenvalue with the 
next-highest magnitude will not be negligible if the difference between the magnitudes of these eigenenergies is not 
sufficiently large. Since the convergence of the power method is generally much poorer than other methods we will 
discuss below, it is generally not used in practice. However, the approach is very simple to implement and is very 
memory efficient because only the two vectors |v„) and |v„_i) must be stored in memory. 
The subspace generated by the sequence of steps in the power method, 

{|vo),//|vo),// 2 |vo),...,i/"|v )} (18) 
is called the n th Krylov space and is the starting point for the other procedures. 



2.3.1. The Lanczos Method 



In the Lanczos method 1 29], the Hamiltonian is projected onto the Krylov subspace using a basis generated by 
orthonormalizing the sequence of vectors Jl 8I > to each other as they are generated. This results in a basis in which the 
matrix representation of the Hamiltonian becomes tridiagonal. The basic algorithm is as follows: 

0) Choose an initial state \uq) which can be represented in the system's many-body basis \n), and which has finite 
overlap with the groundstate of the Hamiltonian. This can be done by taking {|n)} as a vector with random 
entries. 

1) Generate the states of the Lanczos basis using the recursion relation 



\u n +i) 
where a n 

and bi 
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(19) 



with bo = and =0. Note that the Lanczos vectors in this formulation of the recursion are not normalized. 

2) Check if the stopping criterion (w„+i < £ is fulfilled. 
If yes: carry out step 4) and then halt. 

If no: continue. 

3) Repeat starting with 1) until n=M (the maximum dimension). 

4) If the stopping criterion is fulfilled, diagonalize the resulting tridiagonal matrix 
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(20) 



using the QL algorithm. Note that the off-diagonal elements have the value bi and not bj for normalized 
Lanczos vectors |m,). The diagonalization yields eigenvalues Eq, E„, and eigenstates |y/b), ■ ■ ■ , |Vn) which 
are represented in the Lanczos basis. 
5) In order to avoid storing all of the basis states \uq) , \u„), the procedure can be repeated (starting with the same 
initial vector \uo)) to calculate the eigenstates represented in the original many-body basis. This is necessary in 
order to be able to calculate properties dependent on the wavefunction, such as quantum mechanical observables. 
One obtains the coefficients a„ of the eigenstate in the original many -body basis \n) by carrying out the basis 
transformation 



i 

= EE c <' ("I «'■)!") 



(21) 



with a„ = £c ( - (n\ui). 

i 

The algorithm is memory efficient, since only the 3 vectors |m„-i), \u n ), and \u n+ \) need to be stored at once. Note 
that it is also possible to formulate the algorithm in a way such that only 2 vectors need to be stored at the cost of 
complicating the algorithm moderately; the memory efficiency is then the same as the power method. As is typically 
the case in iterative diagonalization procedures, the most time-consuming step is carrying out the multiplication H\u n ), 
which should be implemented as efficiently as possible, either as described in Sec. 12. II or using other sparse-matrix 
multiplication routines. The time needed to perform the other steps of the algorithm is generally negligible when is 
realistically large. Therefore, optimizing the routine performing H\ y/) is crucial to the efficiency of the implementation. 



The Lanczos procedure results in a variational approximation to the extremal eigenvalue which usually attains 
quite high accuracy after a number of iterations much smaller than the dimension of the Hilbert space. Typically, 100 
recursion steps or less are sufficient to attain convergence to almost machine precision for ground state properties 
l3(ill . As is evident because the Lanczos method is based on the power method, convergence to extremal eigenvalues 
occurs first; additional iterations are then necessary to obtain converged excited states. The algorithm is generally 
considered to be a standard method which can be robustly applied to a wide spectrum of systems. Nevertheless, 
there are two technical problems which require some care to be taken. First, the convergence of excited states can be 
irregular; in particular apparent convergence to a particular value can occur for some range of iterative steps, followed 
by a relatively abrupt change to another, substantionally lower value. It is therefore important to carry out sufficient 
iterations for the higher excited states. Generally, the number of iterations required to obtain convergence becomes 
larger for higher excited states. 

The second technical problem is the appearance of so-called "ghost" eigenvalues, i.e., spurious eigenvalues which 
cannot be mapped to eigenvalues of the original Hamiltonian. The origin of such "ghost" eigenvalues can be traced to 
the loss of the orthogonality of the Lanczos vectors \u„) due to finite machine precision. This is an intrinsic limitation 
of the algorithm. For this reason, the much more stable Householder algorithm 1 17] rather than the Lanczos procedure 
is generally used to transform an entire matrix to tridiagonal form. However, because of the good convergence to 
the ground state and the algorithm's memory efficiency, the Lanczos method is nevertheless widely used for the 
investigation of quantum lattice systems with short-ranged coupling. 

With some extra effort, it is also possible to overcome the loss of orthogonality which is a basic drawback of 
the algorithm. The most straightforward solution is to reorthogonalize the Lanczos vectors relative to each other 
using a modified Gram-Schmidt procedure. However, this requires all vectors \u n ) to be stored in memory, so that 
the advantage of memory efficiency is lost. However, an appropriately chosen partial reorthogonalization is also 
sufficient; for details see Ref. H25L \26L 12711 . Cullum and Willoughby |28, 31] developed a method to eliminate ghost 
states without reorthogonalizing the Lanczos vectors. In their approach, the eigenvalues of the resulting tridiagonal 
matrix T„ are compared to the ones of a similar matrix T„, which can be obtained by deleting the first row and column 
of T„. This gives a heuristic criterion for the elimination of spurious eigenvalues: since the ghost eigenvalues are 
generated by roundoff errors, they do not depend on the initial state |vq) and will be the same for both matrices. After 
sufficiently many iterations, ghost eigenvalues will converge towards true eigenvalues of the original matrix H. Thus, 
every multiple eigenvalue of T„ is not a ghost and every unique eigenvalue which is not an eigenvalue of T„ is a true 
eigenvalue of H. This approach is as memory-efficient as the original Lanczos algorithm and approximately as fast, 
but generally yields the wrong multiplicity for the eigenvalues. 

A variant of the algorithm which is easy to implement is the "modified" Lanczos method i30ll . In this approach, 
the recursion procedure is terminated after two steps, i.e., only two Lanczos vectors are considered and the resulting 
2x2 matrix is diagonalized. The resultant eigenvector is taken as the starting point for a new 2x2 Lanzcos procedure; 
this process is then repeated until convergence (i.e., the change in |Ao| is sufficiently small) is achieved. The modified 
Lanczos method has only limited usefulness: the convergence is only marginally better than the power method and 
excited states are rather difficult to obtain. However, the idea of restarting the Lanczos procedure is often used in 
practical implementations, i.e., after ~ 10 to ~ 100 Lanczos iterations, the resulting tridiagonal matrix is diagonalized 
and the extremal eigenstate is used as starting vector for a new Lanczos procedure. Further important variants of the 
algorithm are the implicitly restarted Lanczos method and the Band- or Block-Lanczos method. The latter puts the 
"modified" Lanczos variant on a more systematic footing and deals with the problem that it is not a priori known how 
many iterations are needed for convergence by limiting the number of steps and then restarting the Lanczos-procedure 
with an appropriate, better initial vector obtained by taking into account the outcome of the previous Lanczos run. The 
Block Lanczos approach is useful when degeneracies are present in the desired eigenvalues and all eigenstates and 
eigenvalues of the corresponding subspace must be obtained. Rather than starting with a single initial state vector \ uq), 

a set of p orthonormal state vectors |wq 1 ' ) ) , \ u^), . . . , |hq ) is used to initialize the Lanczos procedure. For an eigenstate 
with degeneracy g c /, p> g c / should be chosen in order to fully resolve the degenerate subspace. The Lanczos procedure 
is then performed on this subspace and a block-tridiagonal matrix is obtained. For details and for other variants of the 
Lanczos procedure, see e.g., Ref. Il27ll . 

The generalization of the Lanczos method to non-hermitian operators is the Arnoldi method, see, e.g., Ref. 12711 . For 
such problems, a Rrylov-space approach similar to the Lanczos procedure is used to reduce a general matrix to upper 
Hessenberg form. 

A number of software packages are available which provide implementations of iterative diagonalization routines, 
such as ARPACK 13211 or the IETL 13311 . a part of the ALPS -library |9]. In order to use this software, a routine which 
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FIGURE 4. Typical convergence behavior of the Lanczos algorithm. In this example, convergence to the ground-state is already 
reached after approximately 10 iterations, while additional iterations are required to reach convergence for the excited states. 
"Ghost" eigenvalues appear and converge to a real eigenvalue as additional iterations are performed, leading to erroneous multi- 
plicity at the completion of the calculation. 

performs the multiplication H\y) must be defined by the user; an effient implementation is crucial to the overall 
efficiency of the algorithm. 

There have been many investigations using the Lanczos procedure. Examples for recent work using ED are 
investigations treating frustrated quantum magnets, Ref. [34], two-leg ladder systems, Ref. 13511 . transport through 
molecules and nanodevices 13611 - and atomic gases on optical lattices 13711 - 



2.3.2. Davidson and Jacobi- Davidson 

The common idea of iterative diagonalization methods is to project the matrix to be diagonalized onto a subspace 
much smaller than the complete Hilbert space. The subspace, spanned by a set of orthonormal states is then 

expanded in a stepwise manner so that the approximation to the extremal eigenstates improves. At particular points in 
or at the end of the procedure, the representation of the Hamiltonian matrix in this subspace is diagonalized. The 
resulting extremal eigenvalue X\ is called the Ritz-value and the corresponding eigenvector the Ritz-vector |t//^). 
According to the Ritz variational principle, is always an upper bound to the real ground state energy, see, e.g., 
Ref. 1 38]. The error in applying H to the eigenvector |%) associated with the Ritz-value /Lj, is approximated by the 
residual vector 

\r k )=H\y k )-l k \%) . (22) 

(This expression would be exact if A* were replaced by the exact eigenvalue X^-) In the Lanczos procedure, the 
recursion is formulated so that the subspace is expanded by the component of the residual vector |r^.) orthogonal 
to the subspace. 

In 1975, Davidson formulated an alternate iterative algorithm in which the subspace is expanded in the following 
way 1 39]. The exact correction to the Ritz-vector is given by 

k) = \¥)-\¥k) 

so that 

(H-Afcl)|z) = -(H-Afcl)|$fc>. (23) 



Thus, solving 

(R-X k l)\z) = -\r k ) (24) 

would lead to the exact correction to This amounts to inverse iteration (see, e.g., Ref. (vfti ). However, the exact 
eigenvalue A# is not known, and the numerical solution of this linear system is of comparable numerical difficulty to 
the entire iterative diagonalization. Davidson's idea was to approximate the correction vector \z) by 

|?) = -(D-A,l)" 1 kA-), (25) 

where the diagonal matrix D contains the diagonal elements of H. This is a good approximation if H is diagonally 
dominant. If D were replaced by the unit matrix 1, the Davidson algorithm would be equivalent to the Block 
Lanczos procedure. Therefore, if the diagonal elements of H were all the same, both methods would have the same 
performance. For many problems, however, this variation of the diagonal elements is important and the Davidson 
algorithm converges more rapidly. 

In its original formulation, the Davidson algorithm follows the procedure [39]: 

0) If the kth eigenvalue is to be obtained, choose a subspace of I > k orthonormal vectors |vi), |v2), . . . , |v/). In the 
following, the matrix B is the matrix containing these vectors as columns. 

1) i) Form and save the vectors H\v\), H\v2), H\v[). In the following, the matrix A is the matrix containing 

these vectors as columns. 

ii) Construct the matrix A = (vi\H\vj) and diagonalize it, obtaining the kth eigenvalue hi and the correspond- 
ing eigenvector |otj ). The upper index (/) denotes that this eigenpair was obtained by keeping I vectors for 
building the matrix A. 

2) Form the residual vector corresponding to the kth eigenvector, 

\ qi ) = (A-Xi l) B)\a { k ' ] ). 

3) Calculate the norm | \qi) |. If | \qi) \ < e, accept this eigenpair, otherwise continue. 

4) Compute the correction vector |w/) = (D — \qi), where D is a matrix containing the diagonal elements of 
A and I is the unit matrix. Orthonormalize against |vi), \v2), ■■■ , \vi) to form |v/ + i). Expand the the matrix 
B by adding |v/+i) as an additional column. 

5) Form and save the vector H|v/ +1 ). Set 1 = 1 + 1 and continue with step 1). 

Although the Davidson algorithm is somewhat more complicated to implement than the Lanczos algorithm, its 
convergence is usually higher order than the Lanczos method and it is more stable. In particular, spurious ghost 
eigenvalues do not appear. The disadvantages compared to Lanczos are that (uj\H\uj) is not tridiagonal, and all the 
\ui) must be kept in order to carry out the explicit orthogonalization in step 4). Similarly to the Lanczos-method, this 
algorithm can fail for particular choices of the initial vector, e.g., if \uq) = |i//o). In practice, one performs a small 
number of Davidson iterations and, if necessary, restarts the procedure using the outcome of the previous iteration as 
initial vectors |vi), . . . , |v/). 

A generalization of the original Davidson approach is given by the Jacobi-Davidson method 1 40] . In this approach, 
one approximates the correction vector \z) by 

(jr-X fc l)|z) = -k)-eh>, (26) 

where £ is chosen so that that (uk\z) = 0. Here the preconditioner J%? is an easily invertible approximation to H. 
The choice of an optimal preconditioner is non-trivial and must be tailored to suit the specific problem treated. The 
Jacobi-Davidson method is more general and flexible than the Davidson method and does not suffer from the lack of 
convergence when \uq) — |yb)- It is (almost) equivalent to the Davidson method when = D. Therefore, it is now 
widely used instead of the Davidson algorithm, especially in quantum chemistry. 

An additional advantage of the Jacobi-Davidson algorithm over the Davidson procedure is that it can be applied to 
generalized eigenvalue problems, 

A\x) =XB\x), (27) 
where A, B are general, complex nxn matrices [40]. 
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name 


notation 


operators 


experiment 


single-particle spectral weight 


A(k,£i)) 


A = c k,a 


photoemission 


structure factor 


Sa(q.«) 


A = S^ 


neutron scattering 


optical conductivity 




A = A 


optics 


4- spin correlation 


R(C0) 


A = I k fi k S k S_ k 


Raman scattering 



Using these iterative diagonalization algorithms, many different quantum many-body systems can be treated. 
However, the methods outlined thus far are primarily useful for calculating properties of the ground-state (highest state) 
and low-lying excited states. Many experimentally interesting properties, such as dynamical correlation functions and 
other frequency-dependent quantities, or finite-temperature properties cannot be directly calculated from the ground 
state and low-lying excited states. Nevertheless, methods based on iterative diagonalization have been developed to 
calculate such properties; these extensions to iterative diagonalization will be discussed in the remainder of this section. 



2.4. Investigation of Dynamical Properties 

Dynamical quantities in general involve excitations at a particular frequency co and are often additionally resolved 
at momentum q. It is clear that extensions to exact diagonalization must generate a set of states that are appropriate 
to describe these energy and momentum scales. This can be done by starting with the ground state obtained by a 
converged iterative diagonalization procedure and applying an appropriate operator or sequence of operators. 

Frequency-dependent quantities can be defined as the Fourier transform of time-dependent correlation functions 

C(t) = -i( WQ \A(t)A\0)\ W o} , (28) 
where the operator A generates the desired correlations. The Fourier transform to frequency space then yields (co > 0) 

C((B + /77) = (v/ |A(a) + /77-// + £ ) _1 A t |v/o) , (29) 
where the inverse operator in the center is termed the resolvent operator. The spectral function then is given by fir! 

I{co) =-- lim ImC(o) + /Tj) . (30) 

K 77^0+ 

Dynamical quantities can be probed in scattering experiments such as photoemission or neutron-scattering. Examples 
of choices of A for some typical spectral functions are displayed in Table|2] 
A spectral function can also be expressed in the Lehmann representation 

/(©)=£ \{Wn\A^\wo)\ 2 8{co-E n +E ), (31) 

n 

which involves a sum over the complete eigenspectrum of the system. Thus, the numerical method must provide 
the poles E„ and the weights (\j/„\A^ \\f/Q) for each eigenstate. While it is possible to calculate the spectral function 
using this formulation in principle, it is generally prohibitive to calculate the entire eigenspectrum of the system; this 
amounts to a complete diagonalization of the Hamiltonian. 

An alternative approach is to apply Eq. (I30i directly. In this case, the resolvent operator must either be calculated 
directly or approximated. 



2.4.1. Krylov Space (or Continued Fraction) Method 

Within the Lanczos method, it is convenient to represent operators in a Krylov basis. The frequency-dependent 
correlation function d29i can be expanded in such a basis if a Lanczos procedure is carried out starting with the initial 
vector 

M = n 1 4 4 1 1 / t|V/0) - (32) 



Within the Lanzcos basis generated by H\hq), H 2 \un), the resolvent operator can be expressed in terms of the 
Lanczos coefficients a„ and b„ as a continued fraction 1 30 ] . The correlation function then has the form 

C(z = co + iri+E ) = -n — ■ (33) 

Therefore, each set of Lanczos coefficients a„, b„ generates an additional pole of the frequency-dependent correlation 
function. In practice, the procedure is truncated after a number of poles appropriate to the size of the system and 
required resolution of the spectral function has been obtained. The weight of each pole typically decreases rapidly 
with number of steps n so that a moderate number of iterations is usually sufficient in order to obtain the desired 
quantities. 

The low energy part of the spectrum can also be obtained by directly using the Lehmann representation of the 
spectral function Eq. J3 11 . Truncation after a finite number of terms yields a finite number of poles. However, a fairly 
large number of eigenstates would have to be calculated directly, which is awkward within iterative diagonalization. 
Additionally, it is not clear that the weight of each subsequent term would dimish sufficiently rapidly to justify a 
truncation after only a small proportion of the eigenstates are obtained. The resulting spectral function using this 
approach is a series of sharp delta peaks located at the poles. In order to obtain a continuous spectrum and to be able 
to compare to experimental data, it is necessary to introduce a broadening of the peaks by calculating the convolution, 
e.g., with a Lorentzian. 

Note that the appearance of ghost-eigenvalues is not a problem in the calculation of the spectral functions, since the 
matrix elements associated with the ghost eigenvalues are negligibly small This turns out to be true when using 
the Lehmann-representation as well as when applying the continued fraction approach. 

2.4.2. Correction Vector Method 

An alternative approach, suggested by Soos and Ramasesha 14311 . is to apply the resolvent operator directly. The 
vectors 

|0o)=A t |v/ o ), |f ) = ((D + H 1 -H + E Q r 1 \<j> ), (34) 

are calculated directly starting with the ground state | y/b) obtained using iterative diagonalization. The spectral function 
is then given by 

7((B) = -Im<to|0i). (35) 
% 

The correction vector |0i) is obtained by solving the linear system 

(£» + rT7-#+£o)|0i) = |0o>- (36) 

One advantage of this approach are that the spectral weight is calculated exactly for a given frequency CO rather than 
approximated as in the continued fraction approach. In addition, it is possible to obtain nonlinear spectral functions by 
computing higher order correction vectors. This scheme is naturally used in conjunction with the Davidson algorithm, 
in which the correction to the residual vector must be calculated at each step, whereas the Krylov subspace is not easily 
available. 

The disadvantage of this approach is that it generally is more expensive in computer time because the system of 
equations, Eq. d36l >. must be solved for each co desired. 



2.5. Real Time Evolution 

The time evolution of a quantum system is governed by the time-dependent Schrodinger equation 

ihj- t \y(t))=H\y(t)). 

Using the Lanczos-vectors, the time evolution through one interval U c i t | y) can be approximated by i44ll45Tl 
\yr(t + dt)) = e- id, l hH '|y/(f)> wV„(f) (0 v J(t) \y{t)) , 



(37) 
(38) 



where V„ is the matrix containing all the Lanczos vectors \uj). One finds that for a high accuracy only a very small 
Rrylov space is needed; n < 20 is sufficient. Therefore, the matrices in Eq. OHJ are very small, and the time-evolution 
can be computed efficiently. For details see the authors contribution to this topic in this volume, Ref. 1 46] . 

Using this approach, it is possible to study a number of different non-equilibrium properties of strongly correlated 
quantum systems. It can also be used to obtain dynamical properties by computing the Fourier-transform of a time- 
dependent correlation function directly. As described in Ref. 1 46], it is possible to construct a time-evolution scheme 
for the DMRG based on this approach. 



2.6. Calculations at Finite Temperature using Exact Diagonalization 

All approaches presented so far have been formulated for zero temperature. The calculation of finite-temperature 
properties is a more demanding task because thermodynamic expectation values are given by sums such as 

(A) = ^{n\Ae-P H \n) 

with the partition function 

Z = £(n|e-^|n) (39) 

where \n) is an orthonormal basis. As in the calculation of dynamical properties, it is prohibitively expensive to sum 
over all \n) due to the large dimension of the basis. This difficulty is circumvented by a stochastic approach ("stochastic 
sampling of Krylov space") developed by Jaklic and Prelovsek |47]. Their approach enables the calculation of 
thermodynamic properties such as the specific heat, the entropy, or the static susceptibilities of a system, as well 
as its static and dynamic correlation functions. It is therefore a useful tool for comparing to experiments. 
The approach is related to the high temperature series expansion, in which expectation values are given by 

(A) = Z-^t^inl^Aln);, (40) 

n k=Q K ■ 

z = EE^<»l**l»>- (4i) 

Expressions such as (n \H k A\n) can be evaluated using the Lanczos vectors and eigenvalues resulting from a Lanczos 
run: 

(n\H k A\n)=Y,(n\Yn(V>nA\n){el M) ) . (42) 

i= l 

However, the number of vectors needed is still too large. To overcome this, in a second step one introduces a stochastic 
sampling over different Krylov spaces, i.e., the Lanczos procedure is performed repeatedly for a variety of random 
initial vectors, and an average is taken over the samples at the end. This finally leads to the expression 



< A > « ^L^LL^ <«W|AW (43) 



where 
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(44) 



Here £ v denotes summation over symmetry sectors of dimension N s , ^L? denotes the average over/? random starting 

vectors I 1 ?!)''), and £„, is over the m th Lanczos propagation of the corresponding random starting vectors |SP« ). 

This approach is useful if convergence with the number of Lanczos propagations M <C N s and the number of random 
samples R -C N s is sufficiently fast. Due to the relation to the high-T expansion, the T — > °° limit is reproduced 
correctly. One obtains high- to medium-temperature properties in the thermodynamic limit. For finite systems, the 



low-temperature limit is reproduced correctly up to the sampling error l47ll . Aichhorn et al. |48] invented an approach 
which can be used to reduce the sampling errors at low temperature. They use the property that a twofold insertion of 
a Lanczos basis leads to smaller fluctuations at low T. To do this, the expectation value is expressed as 

(A) = ^i(n\e-^ 2 Ae-^ 2 \n), (45) 

Z n 

and the projection to the Lanczos basis is performed. 

This method has been used, e.g., for the investigation of the planar t — J model at finite temperature [49]. 



2.7. Discussion: Exact Diagonalization 

As we have seen, the exact numerical treatment of quantum many body systems is hindered by a number of sig- 
nificant obstacles. However, it is possible to formulate conceptually straightforward numerically exact methods which 
allow the treatment of surprisingly large Hamiltonian matrices, especially when symmetries of the Hamiltonian are 
taken into account. Extensions to the basic method also make it possible to calculate quantities which, in their sim- 
plest formulation, require the calculation of the full eigenspectrum of the Hamiltonian, such as dynamical correlation 
functions, the time evolution and finite temperature properties. Nevertheless, maximum system sizes remain strongly 
limited because of the exponential growth of the many-body Hilbert space with system size. In the following sections, 
we discuss additional numerically exact methods which overcome this restriction, especially for one-dimensional sys- 
tems, albeit with the introduction of additional approximations. These methods are capable of treating much larger 
systems, containing up to as much as several thousand sites, allowing a reliable extrapolation to the thermodynamic 
limit to be performed. However, exact diagonalization nevertheless plays the role of a benchmark for other methods 
because it provides numerically exact results in most cases, and is useful for problems for which the other approaches 
fail. 



3. NUMERICAL RENORMALIZATION GROUP 

The numerical renormalization (NRG) was developed by Wilson as a numerical approach to the single-impurity Kondo 
and Anderson problems, which, in conjunction with various analytic methods, provides an accurate, complete solution 
of the problems |6]. The difficulty in these problems lies in the widely varying energy scales that need to be accurately 
described. For example, Kondo screening occurs at an exponentially small energy scale set by the Kondo temperature. 
Perturbation theory, while providing a good description at higher temperatures, breaks down at the Kondo temperature. 

Here we will be concerned primarily with illustrating the general concepts and features of the method, exploring 
the conditions under which it works or does not work, and relating it to other methods; in particular, to exact 
diagonalization and to the DMRG. More complete treatments of the NRG include Refs. |6, 50, 51], on which the 
material presented here is based. 



3.1. Anderson and Kondo problems 

The single-impurity Anderson model (SIAM) was introduced by Anderson in 1961 to describe a spherically 
symmetric strongly correlated impurity in an uncorrected non-magnetic metal . In lattice form, its Hamiltonian 
can be written 

H AI = e d £n d a + + £ (v M c£ ^ +H.c.) + £ ek^.a (46) 

a k,(T k,a 

where the non-degenerate impurity has energy and n d a = d\,d a is its number operator and d' a the corresponding 
creation operator. The Coulomb interaction U is local and takes place only on the impurity, but the hybridization Vkd 
in general allows scattering from all momentum states. The hybridization function 



A(o) = 7r£|V M | 2 S(o)-£ k ), 

k 



(47) 



along with the density of states of the conduction band, p (to) = £ k 8((Q — £k), determine the behavior of the system. 
Here we follow the usual treatment and consider only the orbitally symmetric case, i.e., V^.d = Vkd and £k = £k- 
This corresponds to taking a completely isotropic electron gas, which can be realized only approximately in a real 
solid. We then need to consider only the s-wave states of the conduction electrons, so that Ck,<r can be replaced by 

c k,l=m=0,a = c k,a- 

A closely related model is the Kondo or s-d exchange model, 

H K =J K S d - so + £ £k4,<T c k,<r (48) 

k.a 

with So = /o (7 5'<t/j/o ^ an d localized Wannier state generated by f a = Y*k c k <r m which the electronic impurity is 
replaced by a localized spin- 1/2 described by the spin operator S</. The Kondo model can be derived as a strong 
couplin g (U ^> V) approximation to the SIAM via the Schrieffer- Wolff transformation in the symmetric case {e d = 
— U /2) flUl . leading to the correspondence Jk = 8V 2 /U. 

In order to bring the SIAM and the Kondo model into a form amenable to numerical treatment, the models are 
mapped onto a linear chain model using a Lanczos tridiagonalization procedure. For the SIAM with orbital symmetry, 
the procedure is carried out as follows. The object is to transform the hybridization into a local term. This can be done 
by defining the localized Wannier state |0, a) = /J a \ V/yac) ( | VA/ac) designates the vacuum) with 

ho = ^Y. V l<dC k ,a (49) 

v k 

and the normalization V = (Y,k\Vkd\ 2 ) 1 ^ 2 ■ One then transforms the kinetic energy of the conduction band, H ( = 
Y.k £k c l „Ci a to the new set of operators, taking it from diagonal to tridiagonal form, by constructing a sequence 
of orthogonal states generated by applying H c : 

|1) = l[// c |0)-|0>vW>] 

Ao 

|»+1) = ^-[H c \n)-\n){n\H c \n)-\n-l){n-l\H c \n))} (50) 

(where we have dropped the spin index a for compactness). Note that this is a rather unusual analytic application of 
what is essentially the Lanzcos procedure of Section l2.3.1l The result is that the SIAM is transformed to the form of a 
semi-infinite tight-binding chain 



H AI - e d 4 + Un d n n d L[ + V £ (/ f J a + H.c 

k.a V 



E [ £ " flafn.a + K {fn,af n+ l,c + H c 
n=Q.a 



(51) 



where, in second quantized notation, the operators f„ a create an electron in state \n). The Lanczos coefficients, 
A,, — (n + 1 \H c \n) and e„ = (n\H c \n), depend on the dispersion £<. and the hybridization function A(o) of the original 
SIAM, and, in general, can be determined, in the worst case numerically, during the Lanczos procedure. In general, 
they do not necessarily fall off with n, a property which is crucial for the convergence of the NRG procedure, as we 
will see below. An analogous procedure can be carried out for the Kondo Hamiltonian l !48l >. starting with the local 
Wannier state |0, d) generated by f Q a = Y.k c k a- 

In order to ensure that the tridiagonal chain formulations of the SIAM or the Kondo model lead to a convergent 
NRG procedure, additional approximations must be made. In particular, a logarithmic discretization of the conduction 
band leads to coefficients that fall off exponentially with n. Here, we will discuss the Kondo model for concreteness, 
but extension to the SIAM is straightforward. In the first step, we take the density of states of the conduction band to 
be a constant, D. As long as there are no divergences at the Fermi level, the low-energy physics should be dominated 
by the constant part, as can be justified by expanding the dispersion in a power series about the Fermi vector kp [6]. 
(Generalizations can be made for other densities of states.) The next is a logarithmic discretization of the conduction 
band [-D < < D] in energy, i.e., a division into intervals D + = [A _ '" +1 ), A~"] and D = [-A~", -A _ '" +1 ^] for the 
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FIGURE 5. logarithmic discretization of the conduction band. 

positive and negative parts, respectively (see Fig.|5J- For each interval, we can expand the electron creation operators 
in a Fourier series, defining 

[k]n 

where [k] n is the set of momentum points in the interval A~(" +1 ) < k < A - ", p is an integer, and (£>,, = 27rA" +1 /(A— 1). 
For an interval A~(" +1 ' < k < A~" in the negative range, operators b\ pa can be defined analogously, It can be shown 
that the operators a\ „ a , b\„ c , and their hermitian conjugates obey the usual anticommutation rules for fermions. For 
H , Eq. ( I48> , the localized state can be expressed as 

/,L=(l-A- 1 ) 1/2 EA-" /2 (4,^ + < ; ,,) (53) 

At this point, the approximation is made to neglect all higher terms in the Fourier series, keeping only one electron 
per logarithmic energy interval, i.e., a\ and bl a . Since the impurity only couples directly to the localized state, 
neglecting these states amounts only to neglecting off-diagonal matrix elements, which can be shown to be proportional 
to (1 — A -1 ). Therefore, the approximation becomes valid in the limit A — ► 1; a more detailed analysis of the errors 
can be found in Ref. |6]. 

After applying the tridiagonalization procedure outlined above to the logarithmically discretized version of the 
Hamiltonian H , we arrive at the effective chain Hamiltonian for the Kondo model 

H K =D £ A-/ 2 (/„y n+lia + H.c.) + 2J K £ fl o 3 ail f (54) 

n=0,O O,^ 

where D = D(l + Ar l )/2. The tight-binding part has the same form as that in Eq. (I46> with £„ = and X„ » 

|(1 +A- 1 )A-"/ 2 for large n 

This form of the tight-binding Hamiltonian for the Kondo model and an analogous one for the SIAM are treatable 
numerically with the NRG. Crucial for the convergence is that the tight-binding coupling decays exponentially with the 
position on the lattice. Physically, this discretization was carefully thought out by Wilson to reflect the exponentially 
small energy scales evident in the behavior of perturbation theory for the Kondo problem 161 15011 . Note that it is crucial 
to adjust the discretization parameter appropriately. If A is too close to unity, the NRG will not converge sufficiently 
quickly, if A is chosen to be too large, the error from the logarithmic discretization becomes too large. In practice, one 
chooses a value around 2, but the extrapolation A — > 1 should, in principle, be carried out. 



3.2. Numerical RG for the Kondo Problem 

3.2.1. Renormalization Group Transformation 

We will now outline the NRG procedure as applied to Hamiltonian \5A\ . The goal is to investigate the behavior of 
the system at a given energy scale by treating a finite system of length L with Hamiltonian 

H L = D £ A-' 2 (fl a f n+h<T + H.c.) + 27 £ fl a S ail f (55) 

n=0,a <T,Ji 

(dropping the superscript/subscript "K" for compactness). Due to the exponentially decaying couplings, a particular 
system size will then describe the energy scale set by Dl = D/ A^ 1 ^ 2 . 



The idea of Wilson was to examine the behavior (i.e., the renormalization group flow) of the lower part of the 
appropriately rescaled eigenvalue spectrum of the sequence of Hamiltonians Hl, Hl+ \, ... numerically. It is convenient 
to rescale the Hamiltonians directly, defining Hl = Hl/Dl- One can relate Hl to Hl+i through the recursion relation 

i/ L+1 =A I / 2 // L + ^(/ L + (J / L+1(7 +H.c.) = <M[H L ]. (56) 

a 

This defines the RG transformation. In principle, this transformation is exact (up to the discretization error associated 
with A). However, if one were to treat each subsequent Hl by numerical diagonalization, the memory and work 
needed would increase exponentially because the number of degrees of freedom is multiplied by four at each step. As 
an approximation, Wilson suggested to keep at most a fixed number m of the lowest-lying eigenstates of Hi at each 
step. For the Kondo problem, the error made at each step can be shown to be of order AT 1 1 2 < 1 0]. 

3.2.2. Numerical Procedure 
The NRG method then proceeds as follows: 

1) Diagonalize numerically, finding the m lowest eigenvalues and the corresponding eigenvectors. 

2) Use the undercomplete similarity transformation Ol formed by the m eigenstates obtained in 1) to transform all 
relevant operators on the L-site system to the new basis. For example, Hi = O'HlOl is a diagonal matrix of 
dimension m, but other operators Al = O t AlOl will not, in general, be diagonal. 

3) Form H^ + i from using the recursion relation ( I56> . i.e., by adding a site to the chain and constructing Hl+\ in 
the expanded product basis. 

4) Repeat l)-3), substituting H L +i forH L . 

This procedure is illustrated schematically in Fig. [6] The procedure can be started with the purely local term Ho 
consisting of the impurity coupled to a single site, but in practice, the relevant first step occurs when the dimension of 
is greater than m. In step 1), one quarter of the eigenstates must be found at a general step for the Kondo problem, 
so that typically "complete diagonalization" algorithms as discussed in Sec. 12. 21 are used. In general, the matrices to 
be diagonalized are block diagonal with respect to the conserved quantum numbers of the system, such as number of 
conduction electrons N and the projection of the total spin S z . As described in Sec. 12. II it is important for numerical 
efficiency to separate these blocks. Once this is done, the matrices to be diagonalized are not particularly sparse. In step 
3), matrix elements of operators linking sites L and L + 1 such as fl a fL+\.a must be constructed in the product basis 
\i, sl) = \i) \sl+i ), with the first ket representing the basis of Hl, and the second the basis of the added site (consisting 
of 4 states for the Kondo problem). The expression for the matrix elements of Hl+\ is 

(i, SL \H L+ i\i',s' L ) = A^ 2 (i\i')(s L+l \s' L+l )E^ 

+ (i\fl j) (sL +l \f L+ha \s' L+l ) (57) 

+ (/|/ L ,j/')(^ + ii/ i t +1 ,j^ + i) 

where the E\ are the eigenvalues of Hl, i runs from 1 to m, and N$ L is the number of electrons in state Sl- 

Note that at a particular point in the procedure, the range of eigenvalues of Hl = DlHl which accurately approxi- 
mates the spectrum of the infinite system is limited. In particular, the accuracy breaks down due to the truncation at 
a scale that is a multiple a of Dl which depends on the value of A and the number of states kept in; typically a is of 
order 10 when m is of order 1000. A lower limit is set by the energy scale Dl'. eigenvalues below Dl are approximated 
more accurately at subsequent steps, i.e., for L' larger than L and smaller Djj. Therefore, the non-rescaled eigenvalues 
El in the range Dl <El< OlDl can be accurately calculated at a particular step. 

3.2.3. Renormalization Group Flow and Fixed Points 

In order to understand the behavior of a system within the renormalization group in general, one searches for fixed 
points of the renormalization group transformation 0, defined by 



&[H*] =H* 



(58) 



FIGURE 6. Schematic depiction of the NRG procedure. 
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FIGURE 7. Renormalization group flow diagram for the Kondo model. 



In analytic variants of the renormalization group, the behavior of the Hamiltonian H is generally parameterized by 
a small number of coupling constants. Finding fixed points then amounts to finding stationary points in the flow 
equations governing these coupling constants. In the NRG, identifying fixed points is little more subtle: in practice, 
when the first p rescaled energy levels are independent of L for a particular (appropriately chosen) range of L, one 
identifies a fixed point. Some insight into the physics of the problem is usually necessary to choose appropriate ranges 
of p and L. Once fixed points are identified, the physicical behavior governed by the fixed point is determined by 
the structure of the low-lying eigenstates. Note that such behavior is not guaranteed for a more general tight-binding 
Hamiltonian. 

For the Kondo problem (and the SIAM in appropriate parameter regimes), two fixed points can be clearly identified; 
one associated with the behavior of model (I54> at J = and one with its behavior at J = °°. In both cases, the behavior 
at the fixed point is easy to understand: for 7 = 0, the impurity is uncoupled to the conduction band and the excitation 
energies are those of the non-interacting tight-binding band extending from to L. For J — °°, the impurity forms 
an infinitely tightly bound singlet with site of the tight-binding chain, effectively removing it from the system. The 
excitation energies relative to the ground state are therefore those of a chain extending from 1 to L, i.e., the excitation 
spectrum for Hi(J = °°) is the same as that for Hl_\{J = 0). An additional complication is that the nature of the 
spectrum for Hi(J = 0) depends on whether L is even or odd; the asymptotic values of the scaled excitation energies 
are different for even and odd L are different, even in the limit of large L. This can be taken into account, however, 
by always applying a sequence of two renormalization group transformations, i.e., by replacing & with in d58l > to 
determine the fixed points and then considering the odd L and even L cases separately. In fact, the crossover from the 
/ = fixed point to the J = °° fixed point amounts to a reversal of the behavior for odd and even L because the zeroth 
lattice site is effectively removed from the chain at the strong coupling fixed point. For details on the structure of the 
excitation energies, see Refs. 

Numerically, one finds that for small L and small but finite J, the structure of the excitations is that of the / = 
fixed point, whereas for large L the structure is that of the J = °o fixed point. Stability analysis, which can be carried 
out analytically near the fixed points, shows that the effective Hamiltonian for the weak-coupling fixed point has a 
marginal operator, indicating that it is unstable, while the effective Hamiltonian for the strong-coupling fixed point has 
no relevant operators, indicating that it is stable |6, 53, 50], in agreement with the behavior observed in the NRG. The 
renormalization group flow from the unstable 7 = fixed point to the stable / = °° fixed point is depicted in Fig.0 



3. 2.4. Calculation of Thermodynamic Properties 

Thermodynamic quantities such as the specific heat or the impurity susceptibility can be easily calculated within 
the NRG if the range of validity of the excitation spectrum is taken into account. Generally, one is interested in the 
impurity contribution to the thermodynamic quantities, derived from the impurity free energy F lmv (T) = —IcbTIiiZ/Zc, 
where Z c is the exactly calculable partition function for the noninteracting conduction band. If the entire eigenvalue 
spectrum were known, the partition function would be given by Z(T) = Trexp(— H K /IcbT). However, at a particular 
stage of the NRG, what one can calculate is the partition function for the truncated lattice 

Z L (T) = Txe- n ^ T = Y,e- E »' kBT = £ e -^/^ _ (59) 



Evidently, Zl{T) can only be a good approximation for Z(T) when the temperature for a particular system size Tj, 
is chosen so that fcg7}, *C OlDl, the largest energy scale accurately described by Hl- We previously argued that the 
minimum energy is set by Dl, the error made in substituting Zl for Z in calculating impurity properties has been more 
rigously estimated to be Di/AkgT in Ref. 1 53]. Therefore, the valid temperature range for a given L is set by 



A -i«^«« 
Dl 



(60) 



where a depends on m and A. In practice, kgTi w Dl is a reasonable choice. 
Experimentally interesting quantities are the impurity specific heat 



and the magnetic susceptibility at zero field due to the impurity 



Zimp 



k R T 



Tr(S[) 2 Tr(Sf) 2 e~ B 'l^ T 



(61) 



(62) 



where S\ and are the z-components of the total spin on the L-site chain with and without the impurity spin, 
respectively. The Hamiltonian H c is that of the noninteracting conduction band. The high- and low-temperature limits, 
as well as the leading behavior around these limits have been calculated analytically. These limits can serve as a check 
of the accuracy of the NRG calculations. For a more complete depiction and discussion of results for thermodynamic 
properties, see Refs. I6l l53l.l5(3l . 



3.2.5. Dynamical Observables and Transport Properties 

Dynamical properties, both at zero and at finite temperature, can also be calculated within the NRG procedure. To 
be concrete, we will discuss perhaps the most experimentally interesting quantity, the impurity spectral function at 
zero temperature (for simplicity): 

A(co) = --lmG((0 + iri) , (63) 
% 

where 

G{t) = -i(y, \Td{t)d\0)\xi/ ) (64) 

is the retarded impurity Green function [c.f. Eqs. ( I28M 31H. For finite L, it is convenient to calculate the spectral 
weight within the Lehmann representation 

A L (co) = ^£|<p|4|0)| 2 S(a-£ p + £ ) + \(0\dUp)\ 2 8(co+E p -E ). (65) 
Ll p 

Since the excitations out of the ground state for Hl are well-represented in the energy range Dl < (0 < OlDl, as 
discussed previously, A {(Dl) ~ A(o) when (0 is chosen to be within this range. In practice, a typical choice is (0 = 2(Ol 
where (Ol = kgTL = Dl- One usually uses Eq. ( 165 \ directly to calculate the spectrum, rather than the more sophisticated 
methods such as the Krylov method outlined in Sec. 12.4. ll or the correction vector method outlined in Sec. l2.4.2l because 
all eigenstates within the required range of excitation energy are available within the NRG procedure; it is then easy to 
calculate the poles and matrix elements within this range. Note that the result obtained is a set of positions and weights 
of 5-functions; in order to compare with continuous experimental spectra, they must be broadened. Typical choices 
for a broadening function are Gaussian or Logarithmic Gaussian distributions of width Dl 1 54, 51]. 

One interesting application of impurity problems comes about in the context of the Dynamical Mean Field Theory 
(DMFT), in which a quantum lattice model such as the Hubbard model is treated in the limit of infinite dimensions 
11551 l56l I57II . The problem can be reduced to that of a generalized SIAM interacting with a bath or host. The 
fully frequency-dependent impurity Green function of this generalized SIAM must then be calculated in order to 
iterate a set of self-consistent mean-field equations. Various methods can be used to solve the impurity problem 
including perturbation theory, exact diagonalization, quantum Monte Carlo, dynamical DMRG (DDMRG), and the 



NRG I57l,l58ll . Dynamical properties at finite temperature can be calculated using similar considerations, as long as 
the additional energy scale kgT is taken into account In particular, the procedure outlined above can 

be used as long as frequencies CO > ksT are considered. For CO < IcbT, additional excitations at higher energies than 
Di < co < OLDl become important. Frequencies in this range can then be handled using a smaller L so that <x>l < KgT. 
Since transport properties such as the resistivity p(T) can be formulated in terms of integrals over frequency of 
frequency-dependent dynamical correlation functions |61], these methods can also be used to calculate them. 



3.3. Numerical RG for Quantum Lattice Problems 

There are a number of quantum lattice problems that have Hamiltonians whose structure is formally similar to the 
tight-binding Hamiltonian for the Kondo model (I54> such as the Hubbard model l|8) or the Heisenberg model (II 01 in 
one dimension. While one could consider carrying out a variation of the NRG procedure on these models in order to 
perform an approximate exact diagonalization on a finite lattice, the physically interesting case is the one in which the 
couplings between nearest-neighbor sites are all equal. This amounts to setting A = 1 in the Kondo case, the point at 
which the convergence of the NRG method breaks down completely because the identification of the length L with 
the energy scale is lost. Nevertheless, adaptations of the NRG procedure were applied to small Hubbard chains |62], 
obtaining an error in the ground-state energy of approximately 10% after 4 renormalization group steps. Later work 
on the spin-1 Heisenberg chain obtained an error of approximately 3% in the ground-state energy of the L = 18 chain 
ll63ll . It is important to note that such calculations are variational so that quantities which characterize the long distance 
behavior such as correlation functions which depend on the wave function can have much larger errors than the error 
in the energy. Finally, an adaptation of the NRG procedure to a two-dimensional noninteracting electron gas was used 
to study the Anderson localization problem in two dimensions |64]. The result obtained was that the system undergoes 
a localization-delocalization transition as a function of disorder strength, a result later discovered to be incorrect: there 
is no transition because two is the lower critical dimension I65tl66ll . 



3.4. Numerical RG for a Noninteracting Particle 

More insight into the breakdown of the NRG procedure for one-dimensional lattice problems with non-decaying 
couplings can be gained by applying a variant of the procedure to the problem of a single particle on a tight-binding 
chain, Eq. 0. We consider the Hamiltonian in the formulation 

H = - L f(\t){e+i\ + \t+i){i\)+2£\i){e\, (66) 

1=1 e=i 

where that state |/) represents an orbital localized on site t. In a matrix representation, this Hamiltonian is tridiagonal 
and is equivalent to a discretized second derivative operator, —d 2 /dx 2 . Note that Hamiltonian (I66> does not include a 
nonzero matrix element between sites 1 and L, so that fixed boundary conditions have been applied to the chain, i.e., 
the wave function is required to vanish at the ends. We modify the Wilson NRG procedure slightly to take into account 
that we are treating a simpler, noninteracting system. There are two significant changes: first, we put together two 
equal-sized systems (called "blocks") rather than just adding a site because the dimension of the Hilbert space grows 
more slowly than in an interacting system: linearly rather than exponentially with the length. Second, the mechanics 
of putting two blocks together is simpler in the interacting system. 

In terms of the procedure outlined in Sec. 13.2.21 the diagonalization of the Hamiltonian, step 1), and the transfor- 
mation of the relevant operators using Ol, a matrix whose columns are the m eigenvectors of with the lowest 
eigenvalues, [step 2)] are carried out as before. The procedure is somewhat simplified because of the less complicated 
structure of Hamiltonian J66l >; for example, there are no S z and N quantum numbers which decouple sectors of H. The 
operators to be transformed are the Hamiltonian for a block Hi = 0]HlOl which is diagonal and = oJHlOz,, 
where represents the connection between the blocks. The procedure is conveniently started at L = 1 for which 
Hi = 2 and Ti = — 1 can be represented as 1 x 1 matrices. Matrix representations of a system of size 2L are then 
formed [step 3)] as 



TABLE 3. Lowest energies after 10 blocking transformations 
for the noninteracting single particle on a 1-D chain with fixed 
boundary conditions, keeping up to m = 8 states. 
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FIGURE 8. The lowest eigenstates of two 8-site blocks (solid circles) and a 16-site block (open squares) for the one-dimensional 
tight-binding model with fixed boundary conditions. 
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(68) 



The procedure is then iterated [step 4)], doubling the size of the system at each step. The matrix to be diagonalized, 
H2L, is of size 2m x 2m at a general step, and the truncation is to an m x m matrix, i.e., one half of the Hilbert space is 
discarded. Note that the procedure is exact for the first few steps, as long as m < L. 

As shown in the first two columns of Table |3j the accuracy of the eigenvalues obtained breaks completely after a 
moderate number of steps. The failure of the NRG for this simple problem as well as the reason for the failure was 
pointed out by Wilson in 1986 llr37ll . For a particular L, the eigenfunctions have the form 



v„ L M 



oc sin 



nlti 
Z+I 



n = 1, 



(69) 



because fixed boundary conditions have been applied to the system. In one iteration of the NRG procedure, a linear 
combination of the for n < m are used to form an approximation to \^{t). As depicted in Fig. [8] the boundary 

conditions that are applied to lead to a non-smooth "dip" in the wave function in the middle of the block of size 2L 
when the NRG transformation is carried out. Clearly, this dip can only be removed by forming a linear combination 
of almost all of the yf„(£). The lesson that is learned, then, is that how the boundaries of the blocks are treated in the 
NRG procedure is crucial: a more general treatment is necessary to formulate a numerically accurate real-space NRG 
procedure for short-ranged quantum lattice models. 



4. FROM THE NRG TO THE DENSITY MATRIX RENORMALIZATION GROUP 



4.1. Better Methods for the Noninteracting Particle 

Perhaps the most obvious way of improving the NRG, at least for the single tight-binding particle, is to apply a 
more general set of boundary conditions to the system to be diagonalized. For example, one can apply fixed (vanishing 
wave function in the continuum limit) or free boundary conditions (vanishing first derivative) by forming the matrix 

(jj/?, fixed rj, \ 

jjfixed.b' J ( 7 °) 

when putting two blocks together, where the boundary condition b can be fixed or free 1631 . For example, 

free. fixed / 1 1 



H" cu = I _i 2 J • (71) 

The Hamiltonian H b 2 '[ can then be diagonalized for all 4 combinations of boundary conditions to obtain a more general 
set of basis functions, albeit an overcomplete set which is not orthogonalized. However, if only a small number 
(e.g., ra/4) of eigenstates from each set {b,b'} are kept in the transformation matrix Ol, the columns of Ol can 
then be orthogonalized numerically using the Gram-Schmidt procedure, resulting in an undercomplete orthogonal 
transformation, as before. The matrices H* 6 ' = OJH^'Ol and Tl = oJh^'Ol must then be transformed to prepare 
for the next step. Note that ft|*' is not diagonal. 

As can be seen in Table [3] this "combination of boundary conditions" method works astoundingly well, producing 
the first few eigenvalues to almost machine accuracy for a 2048-site lattice, keeping only m = 8 states |68|. However, 
the crucial question is whether the technique can be generalized to interacting many-particle systems. Unfortunately, 
no good method has been found to extend these "combination of boundary conditions" methods. The reason is that it 
is not clear that an appropriately general set of boundary conditions can be found for a many -particle wave function, 
which is a complicated function of the coordinates of all the particles. 

However, another method developed in Ref. |68], called the superblock method, can be generalized to interacting 
systems. Rather than applying a general set of boundary conditions to a block of size L, this method forms a new basis 
for H2L and T2L based on the idea that these blocks will eventually make up part of a larger system. In order to do 
this, a "superblock" (with periodic boundary conditions) made up of p > 2 blocks is formed and diagonalized. For 
example, 
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(72) 



The transformation O2L is then made up by projecting the m lowest-lying eigenstates of H? L onto the coordinates of 
the first two blocks, and then orthonormalizing its columns. In other words, if m" (with j = 1 , . . . ,4m) is an eigenvector 
of H^, then a nonorthonormalized column vector of O2L is composed of the first 2m elements, j = 1 , . . . , 2m, of u", 

assuming Hl is an m x m matrix. This new basis is used to transform H2L = 02 L H2l02l and T2L = 02 L T2l02l, as 
defined in (|67j and <|68}. 

The superblock method yields good results for the single-particle tight-binding chain, albeit not as accurate as the 
fixed-free variant of the combination of boundary conditions method 1 68]. As one might expect, the method also 
becomes more accurate as the number of blocks making up the superblock, p, is increased. If p can be made arbitrarily 
large, it can be adjusted so that one is simply carrying out an exact diagonalization of a system of the desired size. 
Since the additional blocks in the superblock allow fluctuations at the boundaries of the blocks whose wave functions 
are used in the renormalization group blocking step, this method is more promising to apply to interacting systems: 
such fluctuations should also treat the boundaries of blocks of a many-body system in a more general way. For a 
single-particle system, the projection of the wave function of the superblock onto the subsystem of interest is trivial: 
it is a one-to-one projection because the Hilbert space of the superblock is a direct sum of that of its subsystems. This 



projection is more complicated for a many-body system: one state of the superblock can project onto many states of a 
subsystem because the Hilbert space grows exponentially, as a direct product of that of the subsystems. 

It is also possible to treat the single particle on a tight-binding chain, Hamiltonian d66i . using an algorithm that 
is a direct adaptation of the DMRG for interacting systems |69]. This algorithm is linear in the system size L and 
involves diagonalizing only a 4 x 4 matrix at each DMRG step 1 70] . Since the single-particle algorithm contains the 
essential ingredients of the DMRG algorithm for interacting systems, we believe that understanding the single-particle 
algorithm is useful in gaining insight into the algorithm for interacting systems. We therefore encourage the reader to 
work through the example implementation of this algorithm 116911 and to examine a similar program which is part of 
the ALPS application library £J. 



4.2. Density Matrix Projection for Interacting Systems 

In the following sections, we will discuss how to generalize the superblock method to interacting systems. This will 
lead to the density matrix renormalization group method (DMRG), one of the most efficient numerically exact methods 
for one-dimensional strongly correlated quantum systems. The reader is encouraged to supplement the discussion 
given here with White's original papers ||7|,|8j] and other introductory and review articles on the method lf7ll|72[|73ll. 
In addition, a collection of papers using the DMRG or treating subjects connected to the DMRG can be found in 
Ref. 1741, a survey is presented in Ref. ll75ll . 

As we have seen in Sec. 14. II the crucial step of the superblock method is a projection of the wave function of the 
superblock onto a subsystem, consisting of two identical blocks for the algorithm for the noninteracting particle. Here 
we consider how to carry out such a projection for a many-body wave function in an optimal way. In order to do this, 
we briefly review the general quantum mechanical description of a system divided into two parts. 

The density matrix is the most general description of a quantum mechanical system because, in contrast to a 
description in terms of the wave function, it can be used to describe a system in a mixed state as well as in a pure 
state 1761 17711 . For a general, mixed state, the density matrix of a system is given by 

p=£C a |¥ a )<¥ a |, (73) 
a 

where the coefficients C a are the weights of the states in the mixed ensemble and are normalized so that Tr p = 1 . 
The density matrix for a system in a pure state would have just one term in the sum. Given p, a subsystem A can be 
described by tracing out the degrees of freedom | j) of the rest of the system, yielding the reduced density matrix 

Pa =Tr U) p. (74) 

The eigenstates of p^ form a complete basis for the subsystem A; its eigenvalues w a give the weight of state a in the 
ensemble and carry the information about the entanglement of the subsystem with the rest of the system. The amount 
of entanglement can be quantified by the mutual quantum information entropy 

S(p) = -Tr L/) (plogp) 

= Wee log Wee . (75) 

a 

This feature gives an interrelation between the DMRG and quantum information theory as discussed in more detail 
in Sec. 16.71 If a system in a pure state | y/) is divided into two pails (A and B), | y/) can be expressed in terms of the 
eigenstates of the reduced density matrices of part A, | <p a ), and part B, \% a ) using the Schmidt decomposition 1 78, 79] 

W) = £\/^|<M \Xa) • (76) 
a 

Here the sum is over the nonzero eigenvalues w a of the density matrices of either part, which can be shown to be 
the same. The behavior of the eigenvalues w a with a depends on the wave function considered and on how the two 
subsystems are chosen. If only a small number (i.e., substantially less than the dimension of the smallest of the two 
density matrices) of the w a are nonzero, | y/) can be represented exactly by the sum over the corresponding states. If 
the w a fall off sufficiently rapidly with a, | y/) can be well approximated by truncating the sum in Eq. (I76> to the m 
eigenstates of the density matrix with the largest eigenvalues. This is the case for, e.g., one-dimensional quantum many- 
body systems with a gapped spectrum, i.e., away from a quantum critical point, for which w a falls off exponentially 



(see the discussion of this topic in Ref. |73]). A useful approximation can be achieved if m can be taken to be much 
smaller than the dimension of the eigenbasis of the density matrix. This approximation can be shown to be optimal in 
the sense of a least-squares minimization of the differences between the exact | y/) and the approximate one and 
is equivalent to the singular value decomposition [17]. Finally, if the w a fall off too slowly or not at all, a truncation 
of the Schmidt decomposition becomes a bad approximation for \y/). The worst case is when all the w a are equal, 
which occurs for a maximally entangled state. A representation optimized for multiple states rather than just | can 
be constructed by including additional states in the density matrix, Eq. (I73> with appropriate weights C a - The Schmidt 
decomposition, Eq. (1761 is then no longer applicable and the approximation to a particular state for a given m will 
become less accurate. 

A measure for the error of this approximation is the sum over the weights of the discarded density-matrix eigenstates, 

N 

Ediscard =X! VV «' ( 77 ) 

m 

with N the dimension of the system's Hilbert space. This is called the discarded weight. Thus, using the basis of density 
matrix eigenstates, an "optimal" description for a quantum many-body system can be found (for details, see Ref. 1 69]). 
In the usual case, the calculation is performed in order to obtain the ground state of the system in a particular symmetry 
sector of the Hilbert space, thus only one target state is kept and the density matrix, Eq. ( 1731 . has just one term in the 
sum. Here and in the following, we will therefore discuss the procedure for a calculation with a single target state; the 
generalization to multiple target states is straightforward. 
The basic procedure to carry out the truncation is then: 

1 . Obtain the ground state | y/) of a finite lattice system using an iterative diagonalization procedure (e.g., the Lanczos 
or Davidson algorithms). 

2. Divide the system in two and obtain the m most important eigenstates of the reduced density matrix of one of the 
subsystems. 

3. Transform the system block into the new (approximate) basis with only m states. 

The DMRG is a numerically implemented variational method based on this truncation. In the following, the imple- 
mentation of this method is discussed. For pedagogical example implementations for the noninteracting case, we refer 
the reader to Refs. j^E^]. 



4.3. DMRG Algorithms 

The goal in formulating DMRG algorithms is to embed an NRG-like iterative buildup and truncation of a system, 
termed "system block", in a larger system, the superblock. As described in the previous section, an iterative diagonal- 
ization is carried out on the superblock to obtain the ground state and possibly some other states, and the eigenstates 
of the corresponding reduced density matrix with the largest weights are used to form a new truncated basis for the 
system block. In order to construct DMRG algorithms, two elements of the procedure must still be formulated: how 
the system block is built up, i.e., how degrees of freedom are added, and how the remainder of the superblock, termed 
"environment" or "environment block", is chosen. In the NRG a single site at a time is added to the system, allowing 
a substantial fraction of the energy eigenstates (1/Ng, where Ni is the number of states on site £) to be kept at each 
step. In the DMRG, it is also important to minimize the number of degrees of freedom, especially since the system to 
be diagonalized, the superblock, contains many additional degrees of freedom in the environment block. Therefore, 
a single site is added to the system block at each step in almost all variants of the DMRG algorithm. There is more 
freedom to choose how the environment block is constructed: the algorithms generally fall into two classes, depending 
on how the superblock evolves with iteration. In the infinite system procedure, the size of the superblock increases at 
each step in a fashion reminiscent of the NRG, and in the finite system procedure, the environment block is chosen 
so that the size of the superblock remains constant, allowing an iterative improvement of the wave function or wave 
functions for one particular finite system. 



H 

FIGURE 9. Superblock configuration for the infinite system algorithm. 



4.3.1. Infinite System Algorithm 

In the infinite system algorithm, the environment block is chosen to be a reflection of the system block, (usually) 
including the added site. Therefore, the superblock grows by two sites at each iteration (as opposed to one site in the 
NRG) and the algorithm can be used to scale towards the infinite-system fixed point, or, as we will see, can be used to 
build up an initial approximation to a system of a particular size, which can then be further improved with the finite 
system procedure. Such additional finite-system iteration is necessary for most systems because the infinite-system 
algorithm is not guaranteed to achieve variational convergence with the number of states kept I69ll80ll . 

The infinite system algorithm for the calculation of the ground state of a one-dimensional reflection symmetric 
lattice proceeds as follows. 

1 . Form a superblock containing L sites which is small enough to be exactly diagonalized. 

2. Diagonalize the superblock Hamiltonian //™ pLr numerically, obtaining only the ground state eigenvalue and 
eigenvector y using the Lanczos or Davidson algorithm. 

3. Form the reduced density matrix p„/ for the current system block from y using p„/ = WtjYi'j' wnere 10 is the 
basis of the system block of size 1 + 1, \ j) the basis of the corresponding environment block, and \jfjj = (i\ 
Note that £' = £ = L/2-l. 

4. Diagonalize pg with a dense matrix diagonalization routine to obtain the m eigenvectors with the largest 
eigenvalues. 

5. Construct the Hamiltonian matrix H^+i of the new system block (i.e. the left block A +s) and other operators 
needed in the course of the iteration (e.g. observables). Transform them to the reduced density matrix eigenbasis 
using H" +1 = O^H" +1 Ol, A^ + i = O^A^+iOl, etc., where the columns of Ol contain the m eigenvectors of p u i 
with the highest eigenvalues, and A^+i is an operator in the system block. 

6. Form a superblock of size L + 2 using H(+\, two single sites and Hf +l - 

7. Repeat starting with step 2, substituting f° r #£ Upel - 

The implementation details for the individual steps will be described in Sec. 15. II In order to obtain an efficient program, 
bases and operators should be decomposed according to the symmetries of the Hamiltonian whenever possible. For an 
introduction to the use of symmetries within the DMRG framework, see Ref. Q. 

This procedure has been formulated to obtain only the ground state, but it is easy to extent it to multiple target 
states by constructing additional states during step 2, either by continuing the diagonalization to obtain excited states 
or by applying operators to \ These additional states are then mixed into the reduced density matrix in step 3 with 
appropriate weights, as discussed in Sec. 14.21 

Here a reflection-symmetric one dimensional system has been assumed. It is possible to generalize the infinite 
system algorithm to non-reflection-symmetric lattices by building up the left half and the right half of the system 
alternately. However, such algorithms have not been particularly important because the finite system algorithm, 
discussed below, can treat such cases and can be generalized to systems that are not one-dimensional chains. 

The results of the procedure are the energies and wave functions obtained in step 2. At this point, matrix elements 
of operators within a state and between states can be calculated, provided that the operators have been formed in the 
appropriate basis, i.e., the same basis in which | y) is represented in step 2. 



o o o o • ■ • 



o o o o 




1 



/ l+l 1+2 



1+3 1+4 



L 



o o o o o • 



• o o o 



FIGURE 10. A schematic depiction of a step in the finite-system algorithm. 



4.3.2. Finite System Algorithm 



The infinite system method has the weakness that the wave function targeted at each step is different because the 
lattice size is different. This can lead to poor convergence or complete lack of convergence with m if the wave function 
changes qualitatively between steps. This can occur for states with some incommensuration with the lattice such as, 
for example, fermions with a nonintegral filling or excited states characterized by a particular wavevector. 

Therefore, an algorithm in which the same finite system is treated at each step is very useful. Instead of convergence 
to an infinite-system fixed point with iteration, there is variational convergence to the wave function or set of wave 
functions for a particular finite system. Such an algorithm can be formulated by choosing a block of appropriate size 
from a previous step as the environment block. 

The finite system algorithm for finding the ground state on a one-dimensional lattice proceeds as follows: 

0. Carry out the infinite system algorithm until the superblock reaches size L, storing H( and the operators needed 
to connect the blocks at each step. 

1. Carry out steps 3-5 of the infinite system algorithm to obtain He+i - Store it. (Now £ ^ £'.) 

2. Form a superblock of size L using two single sites and Hf,_ y The superblock configuration is shown in 
Fig.[K)|where£'=L-£-2. 

3. Repeat steps 1-2 until I = L — 3 (i.e. £' — 1). This is the left to right phase of the algorithm. 

4. Carry out steps 3-5 of the infinite system algorithm, reversing the roles of H( and Hfi, i.e. switch directions to 
build up the right block and obtain B§ +1 using the stored H( as the environment. Store Hf, +l . 

5. Form a superblock of size L using H^\, two single sites and B§ +v 

6. Repeat steps 4-5 until £ — 1. This is the right to left phase of the algorithm. 

7. Repeat starting with step 1. 

This procedure is illustrated schematically in Fig. 1101 

One iteration of the outermost loop is usually called a finite-system iteration or finite-system sweep. If the lattice 
is reflection symmetric, the procedure can be shortened by reversing direction at the reflection symmetric point, 
£ = L/2 — 1, i.e., by using the reflection symmetry to interchange the role of the left and right blocks at this point. In 
this formulation, we have assumed that the infinite system algorithm can be carried out to build up the lattice to the 
desired size and to generate an initial set of environment blocks. If this is not the case, the finite system method can 
still be applied if a reasonable approximation is used for the environment block in the first finite-system sweep. The 
simplest such approximation is to use a null environment block, which is equivalent to the Wilson NRG procedure; a 
better one is to use a few exactly treated sites as the environment. As long as this initial procedure does not lead to a 
system block basis that is too bad, convergence is reached after a relatively small number of finite system iterations 
for most systems, typically between 2 and 10 for one-dimensional systems. 

Note that, as in the infinite system algorithm, it is also possible to obtain results for several target states by mixing 
additional states into the density matrix in step 3. 

We show an example of the behavior of the ground-state energy in the course of a DMRG run for the half-filled 
one-dimensional Hubbard model in Fig. ^2 Note the convergence is clearly variational and that there is a significant 
downwards jump in the energy when the direction is changed at the middle of the chain in the finite system algorithm. 
This is due to a qualitative improvement in the representation of the system as the reflected system block can also be 
used as the environment block. 



L=128, t=l, U=4, (n)=l 




FIGURE 11. Ground-state energy Eq obtained by the DMRG algorithm compared with the exact ground state energy from 
the Bethe Ansatz E^ A for the one-dimensional Hubbard model of length L = 128 with U/t = 4 at half-filling. "Position" refers to 
the position at which a site is added. The infinite system algorithm (not shown) followed by six sweeps of the reflection-symmetric 
version of the finite system algorithm are carried out. 



Since its introduction thirteen years ago, the DMRG has become a standard method for obtaining the ground- 
state properties of one-dimensional short-range quantum lattice models. Among the first and most extensive areas of 
application have been spin chains with half-integral and integral spin l8~fll . as well as frustrated spin chains |82]. 
One-dimensional fermionic systems 1 83, 84, 85] have been treated in many variations. Since the amount of work done 
on such models is quite large, a comprehensive survey is beyond the scope of this pedagogical introduction; we refer 
the reader to Refs. I7ll 172. 17311 for further references. 



5. THE DMRG IN DETAIL 
5.1. Programming Details 

In this section, we discuss the details of implementing the DMRG algorithm, paying particular regard to efficiency. 
We follow Ref. |71] in that we take the one-dimensional Heisenberg model with nearest-neighbor exchange term 

S, • S m = S\S\ +l + \ {S+S M + S;S+ +l ) (78) 

as an example. Any DMRG program will contain the following three crucial elements: 

1. The addition of two blocks (usually the system block and a site). 

2. The multiplication // su P er | \jf) . 

3. The basis transformation of relevant operators on a block to the truncated basis of density matrix eigenstates. 
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FIGURE 12. A schematic depiction of the addition of two blocks. 



When two blocks are added, the part of the Hamiltonian internal to each block as well as all terms connecting the 
two blocks must be constructed. For the example, the terms that must be constructed are 



[Hn] if ,jf = [Hi] u , Sj f + Sff [H 2 ] jf + [Sj] ., [Sj +L 



jj 



Note that the matrix representation of each operator appearing in Eq. i!9\ must be available in the current basis; such 
operators must be constructed and transformed in previous steps. 

The multiplication H super \\j/) is the most time-consuming operation in the DMRG algorithm. The following is an 
efficient procedure to perform this operation within the framework of the DMRG; it should be considered to be the 
standard method. Assuming that the system is composed of two parts, the superblock Hamiltonian can be constructed 
using 

a 

where A™ are the matrix representations of the appropriate operators on the left block, while the matrices B a are 
those on the right block. The index a iterates over all pairs of operators that are needed to construct the superblock 
Hamiltonian, as given in Eq. £791 for the example system. The product with the wave function in the appropriate 
superblock configuration then is 

ErHww/^E^W/w/ • (si) 

i'f a i> / 

For each a, this expression is equivalent to the multiplication of three matrices 

a 

which can be carried out in order m 3 operations. 

In the course of the DMRG procedure, all of the above matrices as well as the operators required to calculate desired 
observables must be transformed into the new truncated basis given by the m density-matrix eigenvectors with largest 
weight. If the transformation matrix 0,j- a is composed of m basis vectors uf,, the operator A,^.,/^ is transformed as 

A aa' = 52 () it.uAi,j i'th' ,■.«• > (83) 

U,i',f 

leading to a reduction of the dimension of the matrix representing A from (m\m.2) x [m\m-i) to m x m. 

In order to make the above procedures as efficient as possible, it is necessary to exploit the system's Abelian 
symmetries. The sums over sets of states can be divided up into sectors corresponding to different quantum numbers. 
Once this is done, only the nonzero parts of the matrix representations of the operators which connect particular 
quantum numbers must be stored and the multiplications can be decomposed into sums over multiplications of these 
non-vanishing pieces. Use of non- Abelian symmetries is also possible, but is substantially more difficult I86. l87ll88il . In 
order to treat large systems and make m as large as possible, it is also important to minimize the use of main memory. 
This can be done by storing the matrix representations of operators not currently needed on secondary storage. In 
particular, relevant operators for blocks that are not needed in the superblock configuration at a given step of the finite 
system procedure, but which will be needed in a subsequent step, need not be retained in primary storage. 



5.2. Measurements 



What are termed "measurements" in the DMRG framework are expectation values of operators calculated within 
a state or between states of the superblock, which are obtained in the iterative diagonalization step. The procedure is 
straightforward provided that the necessary operators are available in the appropriate basis. Given a state yfjj of the 
two-block system, the single-site expectation value (y|S|| y) is given by 

(Y\Si\Y)=ZViA S t\*Vi<j- (84) 

The matrix representation is constructed when the site £ is added to the system block, and must be transformed 
at each subsequent step so that it is available in the basis \i). 

For expectation values of operators on two different sites such as the correlation function (\j/\S z ( S; n \ i/), how the 
operators are constructed depends on whether the two sites I and m are on the same or on different blocks. If they are 
located on different blocks, then the expectation value can be formed using 

(xif\S z A\w)= E VtjWSmhfYi'ji, (85) 

where [5^],-// and [S^] if are the individual single-site operators. However, if I and m are on the same block, the 
expression 

(yl-SSftlv) « L tfj[%i<[S l m ]i>i>>Vi»j (86) 

is incorrect within the approximatation for because the sum over i' extends over the current truncated basis rather 
than a complete set of states. The correct way to calculate the two-site expectation value is to use the relation 

{WlS z m W)= J Lwt j [SlS z m ]aWi'p (87) 
u'.j 

where the operator [SjS z m ]ai has been calculated at the appropriate step; this is correct within the variational approxi- 
mation for Xfij. The general rule is, that compound operators internal to a block must be accumulated as the calculation 
proceeds. In this way, almost all correlation functions as well as more complicated equal-time expectation values 
can be calculated. The calculation of dynamical correlation functions and time-dependent quantities have been made 
possible by recent developments and will be discussed in Sec. l6.3l and Sec. 16.51 respectively. 



5.3. Wave Function Transformations 

The most time-consuming part of the DMRG algorithm is the iterative diagonalization of the superblock Hamilto- 
nian. Here we discuss how to optimize this procedure significantly by reducing the number of steps in the iterative 
diagonalization, in some cases by up to an order of magnitude. 

As discussed in Sec. 12. 31 the key operation and most time-consuming part of any iterative diagonalization procedure 
is the multiplication of the Hamiltonian and an arbitrary wave function, i.e., the operation Z/ super y/ in the DMRG. 
Typically, of the order of 40-100 such multiplications are required to reach convergence. Reducing this number would 
thus directly lead to a proportional speedup of the diagonalization. If the Davidson or Lanczos procedure is started with 
a wave vector that is a good approximation to the desired wave vector, much fewer iterations will be required. Since 
an approximation to the same system is treated at each step of the finite system algorithm, an obvious starting point is 
the result | V^o) °f a previous finite system step. However, this wave function is not in an appropriate basis for H super 
because it was obtained using a different superblock configuration. In order to be able to perform the multiplication 
fl| S " p 1 er | Vo), the wave function must be transformed from the basis at step £ to a basis suitable to describe the system 
configuration at step £ + 1 . 

At step i, a state in the superblock basis is given by 

\a e s M s e+2 Pt+3) = \ae)®\st + i)®\s t+ 2)®\Pt + 3) , (88) 

where |cfy) is the basis of the left block containing sites 1,...,£, |s^+i) and \se + 2) are the bases of the single sites, £+1 
and £ + 2, and IPc+i) is the basis of the right block in the 4-block description of the superblock, as depicted in Fig. 
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FIGURE 13. Superblock configuration relevant for the wave function transformation. 



IT3l Assuming the algorithm is building up the system block from left to right, these states must be transformed to the 
configuration of the superblock at step I + 1 with basis 

1^+1^+2^+3^+4) ■ (89) 

This transformation is performed in two steps. The left block is transformed from the original product basis {\(Xe) ® 
|s£+i)} to the effective density matrix basis obtained at the end of the previous DMRG step, {|a^+i)} using 

\a e+l )= £ L e+1 [st + i} ae+1 , at \a e )®\s e+1 ) , (90) 
n+i' a e 

where the transformation matrix L I+l [s/+i]a /+1 ,a f contains the density matrix eigenvectors u"^ a( and is a simple 
rearrangement of the matrix elements of the transformation matrix S(+l a c ;a t+l used in Eq. (I83i for the transformation 
of operators to the new (truncated) density matrix basis. Similarly, for the right basis one defines 

|fc +3 >= £ R e+3 [s M ]p M ,p e+4 W + 3)®\P M ) ■ (91) 

Transformation matrices similar to L l+1 [s/+i]a /+1 ,a f and R t+i [ s e+3]p t+3 ,p t+4 were introduced by Ostlund and Rommer 
in Ref. |89], where they show that the resultant wave function is a matrix product state. Such states will be discussed 
in more detail in Sec. 16.61 

To perform the wave function transformation needed for the left-to-right part of the DMRG, we expand the 
superblock wave function at step I as 

\y)= £ yr(ttt,st+i,Si+2,Pe+3) laese+ise+iPi+i) ■ (92) 

The basis transformation is formally performed by inserting L« f+1 \ &e+i) (c^+i I- Since there is a truncation, this is 
only an approximation, 

£|o%+i)(o%+i|wl. (93) 

at+i 

The coefficients of the wave function in the new basis therefore become 

£ L e+l [s e+1 ] ae+l ,a e V / («^^+i I ^+2 J ft+ 3 )^' +3 [^+ 3 ] ft+3A+4 . 

It is convenient to perform the procedure in two steps: 

1 . Form the intermediate result 

\ir(at + i,s t+ 2,Pe +3 ) = £ L e+l [s M } at+uat y/(ae,s e+u s e+2 ,l5 e+ 3) ■ (94) 

a t ,s t+l 

2. Then form 

^(^+1,^+2,^+3,^+4) = £ \jf(oc e+ us e+ 2,Pe+3)R t+i [se+3]p M ,p M ■ (95) 

Pe+3 



An analogous transformation is used for a step in the right to left sweep. 

The wave function transformations are relatively inexpensive in CPU time and in memory compared to other 
steps of the DMRG procedure. In addition to leading to faster convergence in the iterative diagonalization due to 
the good starting vector, wave function transformations also allow the convergence criterium of the Davidson- or 
Lanczos-algorithm to be relaxed so that the variational error in the iterative diaogonalization is comparable to the 
variational error of the DMRG truncation, saving additional iterations. Normally such a relaxation would lead to 
the strong possibility of convergence to a qualitatively incorrect state. However, since the initial state comes from a 
diagonalization of the same finite lattice, there is little danger of this occuring. 

Implementing this transformation requires saving the transformation matrices L[^] and R[s(] at every finite-system 
step, which was not necessary in the original formulation of the algorithm. In an efficient implementation of the DMRG 
algorithm, the transformation matrices (as well as the matrices needed to describe a block) are stored on hard disk. 
This additionally makes it possible to reconstruct all operators after the final ground state wave function | i/Zq) has been 
obtained, which saves memory during the DMRG run when many measurements are made. 



5.4. Extensions to Higher Dimension 

5.4.7. Real-Space Algorithm in Two Dimensions 

The extension of the DMRG algorithm to higher dimensional systems, i.e., to two or three-dimensional quantum 
systems or to three-dimensional classical systems is a difficult problem. The most straightforward, and, until now, most 
extensively used, way of extending the DMRG algorithm to, for example, systems of coupled chains, is to simply fold 
the one-dimensional algorithm into the two-dimensional lattice, as depicted in Fig.^]| 83, 90]. In effect, one is treating 
a one-dimensional lattice with longer-range interactions generated by the coupling between rows. Note that one site 
is added to the system block at each step as in the one-dimensional algorithm. One could consider instead adding 
port ion of the lattice that is more appropriate for preserving the two-dimensional symmetry such as a row of sites 
19111 : however, this leads to a prohibitive explosion in the number of states in the superblock basis unless additional 
approximations are made. The infinite-system algorithm outlined above cannot be straightforwardly extended to the 
two-dimensional lattice because there is no reflection symmetry between the system and environment blocks at a 
particular point. However, once the system has been built up so that a set of system blocks has been generated, the 
finite-system algorithm can be used with no problems, even taking into account reflection symmetry, if it is present. 
As long as this set of system blocks is a reasonable approximation, the convergence in the number of finite-system 
sweeps is rapid for most systems. 

In order to carry out the buildup and form the first generation of system blocks, different schemes can be used. 
The simplest is to not include the environment at all, which amounts to carrying out the Wilson procedure. However, 
this algorithm has serious drawbacks for quantum lattice problems, as discussed above. In addition, for systems of 
fermionic or bosonic particles at fillings with a non-integral number of particles per site, a chemical potential must be 
introduced and carefully adjusted so that the average particle density is appropriate. An improvement on this scheme is 
to take the environment block to be a small number of exactly treated sites, typically 3 to 8, depending on the system. 
This scheme avoids the problems in the Wilson procedure, is flexible and easy to implement, and generates sufficiently 
good system blocks in most cases. Another possibility is to use a hybrid scheme in which the finite-system algorithm 
is run on a lattice smaller than the target size, and the system size is increased a row at a time when an appropriately 
sized block are available 1 83 ] . (This scheme seems to have been rediscovered in Ref. 19211 ,)_Additional improvements 
to two-dimensional algorithms include adding symmetry-adapted bands rather than sites 19 311 and a scheme in which 
a the choice of particular diagonal path through a two-dimensional lattice makes it possible to formulate both infinite- 
system and finite-system algorithms that add only one site at a time and break up the lattice more symmetrically 
& 

However, all of these schemes suffer from the fundamental limitation that the boundary between the system and 
environment blocks is relatively long (proportional to the linear dimension in two dimensions) with many interaction 
terms along the boundary. In practice, this leads to a large entanglement between the system and environment blocks 
and therefore a slower fall-off of the eigenvalues of the density matrix, requiring many more states to be kept to 
attain a given accuracy. This scaling of the required number of states has been studied systematically for the case 
of two-dimensional noninteracting spinless fermions, where it is found that the number scales exponentially with the 
linear dimension of the lattice II90L 19511 . It is not completely clear that this exponential scaling is fundamental for 
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FIGURE 14. Superblock scheme for a two-dimensional lattice. 



all systems; in fact, noninteracting particles are probably the worst case because of the absence of length scales in 
the wave functions and because of the highly degenerate gapless excitation spectrum. For two-dimensional gapped 
systems, Chung and Peschel |96] have shown that the scale of exponential decay of the density-matrix eigenvalues for 
a system of interacting harmonic oscillators diverges linearly with system width; an argument exists that this behavior 
is generic |93]. In practice, multichain systems up to a particular width that depends on the system studied can be 
treated with sufficient accuracy to obtain well-controlled results. 

The DMRG has been used to study a number of multi-chain and two-dimensional systems, including the multichain 
Heisenberg model |97], the frustrated Heisenberg model in one dimension 1 82] and in tw o dimensions I93| |. the two-leg 
Hubbard ladder l98lTcaVAOQ 0, and the multichain and two-dimensional t-J model fTool fToil f702l 1 1 03ll . 

Recently, a variational scheme related to the DMRG and based on a tensor-based generalization of matrix product 
states has been argued to provide a better 1 104] representation of two-dimensional wave functions. While the prelimary 
calculations look promising, the usefulness of this scheme in practice relative to the existent DMRG algorithms has 
not yet been explored. 



5.4.2. Momentum-Space DMRG 

Another approach to treating higher dimensional systems has been to work with a momentum-space formulation of 
a local Hamiltonian such as the Hubbard model. The Hamiltonian of the Hubbard model in momentum space reads 

H = L £ k c L c ka + Jf E C l- q fk+cil C kl C pl > < % ) 
ko pkq 

where c\ = 1 / y/N e ' kt (: l a is the Fourier-transformed electron creation operator and N is the number of lattice sites. 
DMRG methods were first applied to this Hamiltonian by White in 1993 in unpublished work 1 105] and subsequently 
independently implemented by Xiang with some additional algorithmic improvements 1 106]. 

The basic idea is to carry out the standard finite-system DMRG algorithm (see Sec. I4.3.2> using the single-particle 
momenta as lattice points. The motivation for treating the momentum-space model is that all information about the 
dimensionality and lattice structure in Eq. d96i is contained in the first term which is diagonal. Therefore, the hope is 
to avoid the drastic loss in accuracy with dimension present in the real-space DMRG. The momentum-space basis is 
also a natural starting point for treating low-energy excitations relative to the filled Fermi sea. In addition, unlike in 
the real-space algorithm, it is easy to preserve and take advantage of the translational symmetry of the system, which 
translates to the conservation of total momentum in Hamiltonian ( I96> . This can be used both to reduce the dimension 
of the Hilbert space in the diagonalization step and to calculate momentum-dependent observables. The disadvantage 



of the momentum-space formulation lies in the fact that the interaction term, the second term in Eq. ( I96K is highly 
non-local, which leads to both technical and fundamental problems. 

There are three significant technical problems: the first is that, unlike for a short-range interactions like those in 
the Hubbard model in real space, there are formally a large number of terms in the interaction; in this case N 3 . The 
impact of this problem can be reduced by factorizing the interacting into sums of terms that are internal to one block. 
As shown in Ref. 110611 . this reduces the number of terms linking the parts of a two-part system to 6N. Second, as in 
the two-dimensional algorithm described above, it is not possible to carry out the infinite-system algorithm to build up 
the lattice due to the lack of reflection symmetry at each step. However, a similar problem is present in the simplest 
variant of the real-space DMRG for two dimensions, as discussed in Sec. 15.4. fl some of the same techniques can be 
used here, such as taking a null environment block (equivalent to the Wilson procedure) or taking a small number of 
exactly treated sites as the environment. Third, it is not immediately clear how to order the momentum sites optimally 
in the DMRG algorithm because the interaction term links all sites with the same amplitude. Orde rings in w hich 
momenta that are a similar distance from the Fermi surface are grouped together seem to be favorable 1 107l ll08ll . but 
the issue is still being explored [ 1 09 ] . This problem also occurs in the DMRG treatment of other nonlocal Hamiltonians 
such as those obtained in quantum chemistry; see Sec. 16.41 

The crucial issue in the usability of the momentum-space DMRG is the convergence of the method, especially 
as compared to the real-space DMRG applied to similar systems. In particular, it is clear that the momentum-space 
DMRG becomes exact in the limit of weak interaction U/t because Hamiltonian J96I becomes diagonal 1 106J. In 
contrast, the real-space DMRG becomes exact for the Hubbard model in the atomic limit t — * 0, but not directly 
in the strong-coupling limit U/t — ► °o. Therefore, for a given dimensionality, there is some interaction strength at 
which the momentum-space DMRG becomes more accurate than the real-space DMR G for t he s ame lattice with the 
same boundary conditions. This occurs at approximately J//f = 8fora4x4 lattice 1106, Il07ll . A comparison of 
the convergence for one- and two-dimensional lattices shows that comparable accuracy is attained when the ratio of 
the interaction to the bandwidth U /W is approximately the same, with some additional small loss of accuracy with 
dimension |107|. 



6. RECENT DEVELOPMENTS IN THE DMRG 
6.1. Classical Transfer Matrices 

In contrast to interacting quantum systems, the ground state of a discrete classical model such as the Ising model 
can be trivially obtained: it is the behavior at finite temperature, characterized for instance by continuous phase 
transitions, that is interesting. The DMRG algorithm is powerful technique for obtaining approximate but very accurate 
information about the extremal eigenstates of an operator that operates on many-body Hilbert space composed of the 
direct products of the Hilbert spaces of its component systems. Transfer matrices for a classical systems have such 
properties and can be used to obtain thermodynamic properties. Therefore, a DMRG treatment of suitable transfer 
matrices, called the Transfer Matrix Renormalization Group (TMRG) can be used to calculate the thermodynamic 
properties of a classical system 111 1011 . In addition, there is a well-known correspondence between a of-dimensional 
quantum and a (d + 1 )-dimensional classical system which makes it possible to calculate the thermodynamic properties 
of some quantum system, as will be discussed in Sec. 16.21 

Consider the two-dimensional Ising model with cylinder geometry so that £ ranges from 1 to L with PBC and n 
ranges from \\oN with OBC (see Fig. I15i. The energy between nearest neighbors i and j is —JsjSj where si = ±1 is 
an Ising spin variable. The symmetrical row transfer matrix in the I direction can then be written as 

r^Vls^exp^!^ |nV(^+ik^+i)|exp(|^4) (97) 
where s = (si, . . . ,Sn) designates a vector of Ising spins, K = J / 'k^T is the dimensionless coupling, and 



W(s'iS' n+l \s n s n+ i) =exp 



— {s' n s n +s„s n+ i +s' n+l s' n + s„ + is' n+l ) 



(98) 



is the Boltzmann weight connecting sites around a plaquette of four nearest neighbors. This matrix has dimension 2' 
and is the object that will be treated with DMRG methods. 
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FIGURE 15. Schematic depiction of the finite row transfer matrix T^> for the two-dimensional Ising model. Open boundary 
condition are applied in the N direction and periodic boundary conditions in the £ direction, which is taken to become infinite. 



The thermodynamic properties can be derived from the partition function Z, which can be written in terms of the 
row transfer matrix as 

Z = Trp=Tr(rW) L = ^(l a ) L (99) 
a 

where the X a are the eigenvalues of T( N \ (We assume that T^ N > is diagonalizable and has positive semi-definite 
eigenvalues.) The trick in calculating Z with a method that can accurately approximate the extremal eigenvalues of a 
finite sparse matrix such as the DMRG or iterative diagonalization is the take the limit L — ► °°, keeping N finite, i.e., 
to make the lattice infinite in one direction only. The partition function can then be approximated as Z ~ Af", where 
X\ is the largest eigenvalue of T^ N \ The (unnormalized) density matrix is given by p = X L |Ai)(Ai | where |Ai) is the 
corresponding eigenvector. 

In order to apply the DMRG method, the system, i.e., T^ N \ must be divided into two parts. If this is done at site M, 
can be written 

rW(s'|s) = r L (s^|s L )w(4^ +1 |^M+i)^(s^|s R ) (ioo) 

where 

r L (sJJs L ) = exp 11 W(s' iS ' n+1 \s„s n+l ) (101) 

is the transfer matrix for the left half-row Sl = (si, . . . ,Sm) and Tr (s' r |sr) is defined analogously for the right half -row. 
The equivalent of the reduced density matrix used in the original DMRG is the left density submatrix defined as 

Pl(s' l \sl) = £p(s^Stf|s L s s ) (102) 

(or the right density submatrix Pr, defined analogously). More than one eigenvalue of is importa nt in the L — > °o 
limit, unlike the eigenvalues of p, but the eigenvalues generally decay rapidly, almost exponentially 1 1 1 lUl 12ll when 
the correlation length is finite. The partition function 

m 

Z = 2>, 2 , (103) 

1=1 



where the (of are the eigenvalues of pi, for the partial system is bounded from above by the full partition function Z, 
and the difference becomes small when m becomes large, or when the the eigenvalues cof fall off sufficiently rapidly. 



The DMRG procedure is analogous to that for the quantum system, using iterative diagonalization to obtain the 
largest eigenstate of rather than the ground state of H as for quantum systems. At each step, the local transfer 
matrix for a single plaquette is added to the finite lattice instead of a quantum site. As for the quantum system, 
the system can be built up to a given row length using the infinite-system algorithm, and then finite-system sweeps can 
be carried out until convergence is reached. One major difference in terms of efficiency compared with the calculation 
for the quantum algorithm is that there are no conserved quantum numbers, reducing the maximum number of states 
m that can be kept for a given amount of numerical effort. 

While all thermodynamic quantities can, in principle, be calculated by taking numerical derivatives of the free 
energy per site / = —ksTlnZ/N « fcgrinAi (where X\ is the variational approximation to Ai obtained from the 
DMRG), second derivatives needed to obtain the susceptibility or specific heat are numerically unstable and it is better 
to use the first derivative of the energy per site, which can be calculated using matrix elements such as (X\\siSj\X), 
which, of course, also yield the short-range spin correlations. The correlation length £, can be calculated from the 
short-range spin correlations if they are small, or directly using £, = l/lnRe^/Ai), where X2 is the second -largest 
eigenvalue. 

Applications of the TMRG for classical row transfer m atrices include the two-dimensional Ising model 1 1 101], the 
spin-3/2 Ising mode l 111 3ll . wetting phenomena | l 14lTl 15ll . the g-state Potts model il 16ll . t he q- state Potts model with 
random exchange | lT7|,Jhe 19-vertex model 11 18ll . and the three-state chiral clock model ill 19ll . For a more extensive 
recent survey, see Ref . 17311 . 

It is also possible to treat Baxter's corner transfer matrices 1 1201 111 ill using the DMRG, leading to the so-called 
Corner Transfer Matrix Renormalization Group (CTMRG) II 2111 . Since corner transfer matrices break up a two- 
dimensional square lattice into quadrants, the partition function is given by the fourth power of the corner transfer 
matrix 

z( w - 1 )=Trp c «Tr(c^) - £ a v 4 , (104) 

' v=l 

and converges variationally to the exact partition function in the thermodynamic limit and as the number of states kept 
m becomes sufficiently large. Here the corner transfer matrix is defined as 

{ s }(ijki} 

and a v are its (largest) eigenvalues. Baxter calculated the free energy per site in the N — > °° limit varia tiona lly 11 2(1 . 
Okunishi pointed out that Baxter's treatment is essentially the infinite-system method of the DMRG Il22ll : Nishino 
and Okunishi formulated a numerical DMRG algorithm based on these ideas, the CTMRG ill 2 ill . It has the advantage 
that the eigenvalues of can be found without diagonalizing a large sparse matrix . The CTMRG has been used to 
find the critical exponents of the two-dimensional Ising model to high precision 112 ill . to investigate the spin-3/2 Ising 
model Il23ll and a 7-configuration vertex model Il24ll . and to study models of two-dimensional self -avoi ding walks 
Efforts are underway to generalize CTMRG-like algorithms to three-dimensional lattices 



6.2. Finite Temperature 

6.2.1. Low Temperature Method 

In the original formulation of the DMRG 1 8], it was clear that a system in a mixed state could the treated within the 
DMRG simply by obtaining excited states of the superblock in addition to the ground state and by including them in 
the system block density matrix: 

a j 

where the weight P a > 0, La^a = 1, and the \y/ a ) are target states. For a system at finite temperatu re, th e natural 
choice of the weight P a might seem to be the Boltzmann weight P a = e~P Ea . Moukouri and Caron Il28ll pointed 
out, however, that P a is only used to form the basis of the system block and does not enter into the calculation of 
expectation values directly. They found that it is more efficient to target a moderate number M (typically 10 to 30) 
of superblock states and weight them equally, i.e., choose P a = 1 /M. In order to calculate thermodynamic properties, 
they form the superblock configuration shown in Fig.ll6lin which the number of states kept in each block is sufficiently 
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FIGURE 16. Superblock configuration consisting of two truncated blocks. 



small so that a complete diagonalization can be carried out. The thermodynamical expectation value of an operator A 
is then approximated by the Boltzmann sum 

{A)=^e-^{ W l B \AWl B ), (107) 
r 

where the set of eigenstates Ey is approximate and incomplete. The assumption is that the truncated basis of Hg 
gives an adequate description of the states in the Boltzmann sum even if they are not included as target states. 
This approximation is clearly better in the low-temperature regime and most likely breaks down at moderately high 
temperature. However, it has the advantage that it is a straightforward extension of the ground-state DMRG and 
therefore can treat the same systems, albeit with a somewhat higher computational cost. This method has been applied 
to S = 1 /2 and S = 3/2 Heisenberg chains Jl2l and to an S = 1 model for Y 2 BaNiQ 5 fl29HT30ll . 



6.2.2. Transfer Matrix Method 



The correspondence between two-dimensional classical systems and one-dimensional quan tum systems a t finit e 
temperature has led to the application of TMRG techniques to quantum transfer matrices il3lL Il32l 
In particular, taking the transfer matrix in the spatial direction leads to a method in which the system is at finite 
temperature, but in the thermodynamic limit. In order to formulate the quantum system as a two-dimensional classical 
lattice model, the partition function is discretized in the imaginary time or temperature direction using a Trotter-Suzuki 
(checkerboard) decomposition 

Z = Tr e-P H = Tr( e - AT/ W- A ^ven)M/2 + ^ (10g) 

where At = j3 /2M is the imaginary time step, 

L-l L-l 

Hodd = 52 n 2t-i,2l and H &ven = £ %,2M-i 
1=1 1=1 

with being the Hamiltonian terms linking site i with site £+1. (We assume that H = H M + H even contains 

only nearest-neighbor interactions.) The error is due to neglecting the commutator in factorizing the exponential. At 
this point a complete set of states \aj) is inserted at every point in space I and imaginary time j so that the partition 
function can be written as a discretized path integral 
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(<% °n + l \e- ATh ^ l \c4t l > ■ (109) 

The resulting lattice in space and imaginary time is shown schematically in Fig.[n] Note that the finite lattice treated 
with the TMRG is in the imaginary-time direction where the boundary conditions are periodic. By defining local 
transfer matrices on a plaquette as 

T(aa\aa) = (\e- ATht -M\) , (110) 
the transfer matrices in the spatial direction for odd and even sites can be written as 
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FIGURE 17. Checkerboard decomposition of the partition function. The quantum transfer matrix in the spacial direction is 
indicated by the region of differently marked squares. 



and 

^even = 12,314,5 ' ' ' T 2M,1 , 

and the partition function for a finite discretization as 

= lim Tr(r odd r e ven) i/2 ■ (111) 

As in the row-transfer TMRG for classical systems, the free energy per site is related to the largest eigenvalue of the 
transfer matrix, k\, 

f M = -^flnX u (112) 

DMRG methods are used to obtain X\. The procedure is a direct adaptation of that for the classical case, except that: 

(i) The detailed expressions for the transfer matrix are complicated by the checkerboard decomposition - the local 
transfer matrices only operate on the filled squares of the checkerboard, and the transfer matrix repeats itself after 
two rows rather than every row. 

(ii) Since the total transfer matrix is given by 7^ M ' = T ^T even and since the commutator [T ^, T even ] ^ in general, 
ji M ) j s no t necessarily symmetric. However, the largest eigenvalue is real for physical cases, and a stable iterative 
diagonalization is possible, for example, with the Arnoldi or unsymmetrical Lanczos methods (see e.g. Ref. 1 27]) 
with additional re(bi)orthogonalization. The left and right eigenvectors and \ \jfg) can thus be calculated and 
are most conveniently chosen to be biorthogonal Vr) = 0- 

(iii) Conserved quantum numbers in the system Hamiltonian lead to corresponding conserved quantities in the transfer 
matrix, which can be used to reduce the dimension of the space treated in the iterative diagonalization. 

(iv) In order to perform basis reduction, the reduced density matrix p = Tr e | y//j) (Vl|> where the trace is over the basis 
of the environment block, is used. This in general unsymmetric matrix must be carefully diagonalized, taking into 
account numerical problems such as the appearance of spurious complex conjugate pairs of eigenvalues and the 
loss of biorthogonality 1 1 3211 - 

As in the classical case, the infinite-system algorithm followed by a number of finite-system sweeps can be carried 
out until convergence is attained. Since we are treating the infinite-system limit of the quantum transfer matrix in 
the spatial direction, the lattice size corresponds to the number of imaginary time steps taken, and the temperature 



P = 1/ksT = MAt. Therefore, it takes more numerical effort to reach a lower temperature and the method should be 
numerically exact at high temperature. Since there is a systematic error associated with the finite imaginary time step 
At, the extrapolation At — » should be carried out to obtain (numerically) exact results. Thermodynamic quantities 
can be calculated as numerical derivatives of the free energy, or can be formed a calculation of the local energy, 
as in the classical case. Dynamical quantities for quantum systems, as discussed in Sec. 12.41 can be calculated at 
finite temperature with two caveats. First, only spatially local or at most short-range expectation values can be easily 
calculated because the transfer matrix only contains spatial information at a range of two sites. Second, the dynamical 
information in the imaginary time direction can be obtained directly from the calculation of the appropriate correlation 
functions on the imaginary-time lattice. However, the experimentally relevant dependence is on real frequency. The 
analytic continuation of imaginary-time data to read frequencies occurs via an inverse Laplace transform, which is an 
ill-posed problem. This issue of how to perform such an analytic continuation numerically and what the limitations are 
has been extensively investigated in the context of quantum Monte Carlo calculations; the Maximum Entropy method 
is the most used technique Il35ll . 

The qua ntum TMRG has been applied to calculate the thermodynamic properties of a number of sys tems ? including 
spin chins 1 136, ll32lll37ll . spin chains with frust ratio n 11381113911140 1. the one-dimensional t-J model 1141»142lll43ll . 
the one-dimensional Kondo lattice model 1 14411 145ll . and a spin-orbit model ll46ll . Dynamic al quantities have been 
calculated for anisotropic spin chains Il47ll . the one-dimensional Kondo model lll48l Il49ll . and the nuclear spin 
relaxation rate in two-leg spin ladders 1 15011 . 



6.3. Dynamics 

As discussed in Sec. 12.41 the calculation of dynamical quantities is a demanding task. In order to obtain the 
dynamical correlation function 

G(k,co) = (\ l fo\Al{co + ir 1 -H)- l A k \xtf ) , (113) 

a ground state calculation is not sufficient. It is necessary to obtain additional states, which in the DMRG must be 
included in the density matrix as target states. The methods vary in which states are targeted, but the state y/o) must 
be targeted in addition to the ground state in all methods. 

The first approach to the calculation of dynamical properties using the DMRG that was developed is the Lanczos 
vector method 1 151]. In this approach, the sequence of Lanczos vectors which span the Rrylov subspace 

IVo), A k |Yb), HA k \Y ), H 2 A k \y ), ... (114) 

are targeted. In this way, the Lanczos coefficients can be obtained and the spectral function can be calculated using 
the continued-fraction representation, Eq. J33l >. This provides the spectral function for the full range of co in one run. 
Usually, one restricts the target vectors to approximately the first 10 vectors of the Rrylov subspace; otherwise, the 
DMRG run would become too expensive because the number of states kept would have to be increased in order to 
obtain the low-lying states sufficiently accurately. However, this number of states is generally not sufficient to reach 
convergence in the continued fraction expansion. In order to overcome this problem, the Lanczos recursion is continued 
using vectors that are not targeted in the density matrix. However, the additional Lanczos vectors are represented in 
a basis of density-matrix eigenvectors which have not been adapted for these vectors. This could lead to uncontrolled 
errors. 

An alternative approach which overcomes this problem is the correction vector method I1521I . It is based on the 
correction vector method f or iter ative diagonalization J43ll discussed in Sec. 12.4.21 and was further developed and 
successfully applied in Ref. lll53ll . In this approach, the target vectors 

|yoM k |i//b), (co + iTj-Hy^Yo) (115) 

are included in the density matrix. The method is more accurate but more costly than the Lanczos vector method, since 
for each value of CO a separate run must be performed. 

The direct inversion of the resolvent operator can be carried out with a procedure such as the conjugate gradient 
method, but this procedure is computationally costly and can be unstable. (These problems can be minimized by a 
better choice of methods for the inversion 1115311 .') A variant of this method was introduced by Jeckelmann 1 154] which 
leads to a more efficient and stable method. The basic observation is that the correction vector minimizes the functional 



W A , r) (co,\i/) = (w\(E + co-H) 2 + ri 2 )\\i/) +r](\j/ \A\\j/} + r](\j/\A\\i/ ) 



(116) 



The spectral weight is then given by 

Wk^OWmin) = -7TT7lmG(k,to) . (117) 

In the calculation, the correction vector is split up into imaginary and real parts 

{(o + iri-H)- l A k \Vo) = \XA{(0 + iri))+i\Y A {(o + iri)) . (118) 

The target-states used in this approach are 

|y/ ),A k |y/ ), \X A (a) + i7])), \Y A (eo + jtj)) . (119) 

The DMRG procedure can be formulated s o that it s imultaneously minimizes the ground-state energy and provides 
the minimum of the functional Wa.ti ((0, w) Il54lll55ll . 

Using this approach, several investigations have been performed. In the following we will outline the procedure by 
outlining the calculation of the single-particle spectral weight for the Hubbard chain 1 156]. The ID Hubbard model 
with open boundary conditions has the Hamiltonian 

L-l L 

H = -t £ (4,a c e+i,a+ c e+i,a c e,c)+ u ll n ^ n U- < 12 °) 

£=l,o v ' i=\ 

The single-particle spectral weight for holes (as obtained in photoemission experiments) is defined by 

A(k,co) = Ilm^ol c\ a H+m ' EQ + ill c k , a \Wo) (121) 

where 

%c = ^Le U Ce,c- d 22 ) 

This definition of c k is problematic because the Fourier transform used is only defined for periodic boundary 
conditions, while a system with open boundary conditions is generally treated by the DMRG for reasons of efficiency 
(see Sec. l4.2> . The solution is to carry out the transform using "particle-in-a-box" eigenstates, 

c ^ = VZTT? sin(H)c ^- (123) 

A comparison of the results obtained with the analytic and exact Bethe-Ansatz method shows excellent agreement. 
Thus, it is possible to obtain momentum-dependent dynamical quantities using open boundary conditions, for which 
the DMRG works best. For this model, relatively large system sizes can be reached (about 200 sites for a Hubbard 
chain on present-day desktop PCs). 



6.4. Quantum Chemistry 

Quantum chemical or ab initio treatments of a solid can be cast into the form of a lattice Hamiltonian once a suitable 
single-particle basis of N orbitals has been chosen. The generic form of the Hamiltonian is then 

H = E Wlcjo + \ E G m E c l c )o> C ko> C ,a- ( 124 ) 
i,j,a ij,k,l a,a' 

where cj a creates an electron in molecular orbital i, is the single -electron integral of molecular orbitals i and 
j, and Gijki is the two-electron integral (Coulomb repulsion) Il57ll . This Hamiltonian can be treated using the 
iterative diagonalization methods described in Sec. 12.31 in quantum chemistry nomenclature, this is termed a "full 
Configuration-Interaction" (full-CI) calculation. Since the Hamiltonian J124l > has a form that is similar to but more 
general than the Hubbard model in momentum-space ( I96> . DMRG methods that are similar to those used in the 



momentum-space DMRG can be applied and similar technical problems crop up. In particular, the molecular orbitals 
take the role of the single-particle momenta with the difference that there is no momentum conservation. A one- 
dimensional path is chosen in this space of molecular orbitals and similar problems occur in the initial buildup of 
the lattice. Both the single-electron term and the two-electron term are non-diagonal in general. The two-electron 
term can be factorized into sums of terms internal to a block as in the momentum-space algorithm; however, the 
minimum number of terms scales as A^ 2 rather than N because of the absence of momentum conservation 1 157]. The 
computational cost of the algorithm scales as (N 4 m 2 + Af 3 m 3 ), where the first term comes from the formation of the 
necessary operators as a site is added to a block and the second is due to multiplication of the 0{N 2 ) representations of 
operators on a block, which comes about either when performing a change of basis on the representation of an op erator , 
or performing the multiplications necessary to form the product H\\ff) needed for the iterative diagonalization Il57tl . 
(Each operation is repeated G(N) times in the finite-system algorithm.) 

As in the momentum-space algorithm, which ordering of the sites, i.e., the molecular orbitals, is optimal is not clear. 
It seems that the ordering can strongly effect both the perfor mance (i.e. , how many finite-system steps are required 
to reach convergence) and the ultimate accuracy of the result 1 158, In particular, some orderings can lead to a 
convergence to higher-lying metastable states. Unfortunately, a general method to find a globally optimal ordering has 
not yet been found. The simplest possibility is to use information from the Hartree-Fock energy and occupation of the 
orbit als to choose th e ordering, typically keeping partially occupied orbitals near each other in the center of the chain 
Il57lll6dll6llfl62il . Another approach is to use information from the single-particle matrix elements ty. In Ref. 
the symmetric reverse Cuthill-Mckee reordering, which attempts to make tu as band-diagonal as possible, was shown 
to lead to faster convergence in m and the number of finite-system sweeps as compared to Hartree-Fock ordering. 
Improvement was also obtained by optimizing the ordering using the Fock term which contains contributions from 
both tij and parts of the interaction Vim I164II . Ideas from Quantum Information Theory also seem to be important: 
Legeza and Solyom 1 109] have shown that the profile along the chain of the entanglement entropy between single sites 
and the rest of the system is important. In particular, orderings which lead to fast convergence tend to have sites with 
high entanglement entropy grouped together near the center of the chain. 

Another possibility to improve the quantum chemical algorithm is to change the choice of single-particle orbitals 
either by choosing a different set of single-particle orbitals at the outset or by performing a canonical transformation 
on the single-particle basis 

j 

where Uu is unitary and and 0, represent single-particle bases. Since such a transformation can change the range 
and topology of connections of both Uj and Vy&/, the performance of the DMRG algorithm can be strongly affected. 
Daul et al. 1 160] applied a procedure in which orbitals unoccupied at the Hartree-Fock level were localized, but did 
not find a significant improvement in the convergence as compared to standard Hartree-Fock orbitals for the HHeH 
molecule. White 1 165] applied a continuous sequence of canonical transformations related to the flow equation method 
Il6dll6ll in order to reduce non-diagonal matrix elements. Calculations for the stretched water molecule as compared 
to full-CI calculations show promising results lll65ll : this method remains to be further explored. 

Thus far, the DMRG for quantum chemistry has been primarily applied to relatively small molecules in order to test 
and develop the method. Full-CI results carried out with the same basis set as the DMRG calculations have provided 
a useful benchmark for the DMRG method. In particular, the water molecule has been treated in the DZP basis with 
8 electrons on 25 molecular orbitals which is small enough so that full-CI results are available 111571 1 60. 163, 108]. 
Agreement with full-CI to within 0.004 mH (Hartree) is possible, keeping up to m = 900 states lll63lll08ll . A TZ2P 
basis with 10 electrons on 41 orbitals, which is too large to treat with full-CI methods has also been treated with the 
DMRG, yielding an estimated accuracy of 0.005m H in t he ground-state energy for m = 5000, approximately a factor 
of three lower than the best coupled cluster method 1 168]. The energy dissociation curves for N? 1 161 ] and H2O 1 162], 
and excited states of the HHeH molecule [ 160] and CH2 [ 108] have been calculated with good accuracy. For LiF, 
the ground and excited states and the dipole moment were calculated to relative accuracies of 10~ 9 , 10~ 7 , and 10" 3 , 
respectively, as compared to full-CI results near the avoided crossing associated with the ionic-neutral crossover of the 
bond 1 164]. 



6.5. Time Evolution 



As discussed in Sec. 12.51 the time evolution of a state in quantum mechanics is governed by the time-dependent 
Schrodinger-equation, Eq. Q, with the formal solution 

\y(t)) = e- m \y(Q)). (125) 

Typically, one is interested in the behavior of the system after a perturbation 

H=H Q +H l ©(f) 

is switched on at t = or after an external operator 

A f |y/ ). 

is applied. In order to calculate the time evolution of the system, one either has to integrate Eq. (0 directly or one 
has to find a suitable approximation for the time-evolution operator exp (—i Ht). I n the context of the DMRG, both 
approaches have been utilized. Recent developments are recapitulated in Ref. 1 169]. 

In the direct integration approach, a Runge-Kutta integration 1 17] of the time evolution was carried out lfT70ll. As 
pointed out in Ref. 1 171], the main difficulty in calculating the time evolution using the DMRG is that the effective 
basis determined at the beginning of the time evolution is not able, in general, to represent the state well at later times 
because it covers a subspace of the system's tot al Hil bert space which is not appropriate to properly represent the state 
at the next time step. As discussed in Refs. Bfl7lll. it is possible that the representation of the time-dependent wave 
functio n very so on becomes quite bad. It is therefore necessary either to mix all time steps |y/(f,)) into the density- 
matrix Il7llll72ll . or to adapt the density matrix during propagation. An approach for adaptiv e time evo lution basing 
on the Trotter-Suzuki decomposition of the time-evolution operator was developed in Refs. Cil[i7i[i7l. In this 
approach, the time-evolution operator is divided into two-site terms which can be applied to the two sites included 
exactly in a typical superblock-configuration. However, the method is restricted to systems with nearest-neighbor 
terms in the Hamiltonian only. A more general scheme to carry out the evolution through one time step that is closely 
related to exact techniques and is combined with an adaption scheme recently proposed by White 1 176] is presented in 
a second contribution of the authors in this volume, Ref. 1 46 ] . In this approach, the time-evolution operator is expanded 
in a Lanczos basis in a similar manner to the time evolution appro ach u sed in the exact diagonalization scheme, as 
discussed in Sec. 12. 51 By using the multi-target method discussed in 1 176, 46], it is possible to efficiently calculate the 
time evolution of systems which can be treated with the DMRG, even if the Hamiltonian contains terms connecting 
sites that are not nearest neighbors, e.g., on ladder systems. However, the how to carry out the targeting optimally is 
still under investigation and will be presented elsewhere 

Since it is possible to treat fairly large systems with the time-dependent DMRG, it is well-suited to investigate 
strongly correlated quantum systems out of equilibrium; such a tool was missing until now. The application of the 
DMRG to such systems opens up new possible applications to a variety of systems in the near future. For recent 
investigations using this method, see Refs. CHS 

and references therein. Note that it is also possible to use the time- 
dependent DMRG to investigate spectral quantities by directly calculating the Fourier transform of time-dependent 
correlation functions. This method should be faster than the DDMRG presented in Sec. l6.3l when a complete spectrum 
is required because information for a wide range of frequencies can be obtained in one DMRG run. However, it is 
likely that the DDMRG will still provide more accurate results for a given frequency 0). 



6.6. Matrix Product States 

In the DMRG procedure, the state of the system block is changed in two steps: First, degrees of freedom are added 
to the system by adding a site (or other subsystem) to the system block, increasing the size of the basis 

\a)t®\st)^\Pl +1 ). (126) 

Second, this expanded system is transformed to a new truncated basis (i.e., the eigenbasis of the density matrix with 
the largest eigenvalues). 

|j3 m )=Z|$ +1 ) (127) 



In fact, both steps can be combined into one transformation, as was done in the wave function transformation discussed 
in Sec. |53] 

|ft+l) = L Ariose) ®|^). (128) 
s(,a 

For each state of the added site s, the matrix Ag a maps from the basis Off describing a system block with I sites to 
Pl+l, the truncated basis for a system block with £+1 sites; typically these bases are the same size. Since the other 
half of the system, i.e., the environment, is built up in a similar fashion (either symmetrically in the infinite-system 
algorithm or in a previous finite-system sweep), its basis can be described as a similar series of transformations. The 
DMRG wave function can then, in general, be written as a sum of such products, leading to the state 18911 

IVmp) = £Tr{A s >A s 2A s 3...A^}|s 1 ,s 2 ,..., iL ) . (129) 

{*<} 

Here the A Sf are m x m matrices, except for A i and Ai which are m-element vectors. (We assume that the same number 
m states are retained at each step; generalization to a variable number of states is straightforward.) The trace is carried 
out over the set of coordinates of the added sites. Generally, a normalization 

^A^(A s <) t = l (130) 

M 

is chosen for the A-matrices. This is a special case of a matri x product state. The DMRG can be viewed as a particular 
variational calculation in the space of such states I89l ll78ll . Note that a matrix product state for periodic boundary 
conditions can be formed if the vectors A\ and Al are replaced by matrices. A translationally invariant state would 
then have all A st = A s equal Il78ll . 

It is useful to consider matrix product states and their properties in more generality. In particular, any state can 
be described as a matrix product state if the dimension of the A* 1 at each step is large enough. However, this is not 
useful in itself: the required dimension can, in principle, grow exponentially with the system size. Nonetheless, some 
states are very compactly described: the Neel state | Tit IT • • •) can t> e trivially described using lxl matrices and the 
Bell state | TTT • • •) ± I III ■ ■ ■) using simple 2x2 matrices. The valence-bond solid state for the spin-1 chain, which 
provides a qualitatively useful description of the ground state fo r the Heis enberg Hamiltonian and is the exact ground 
state of the Affleck-Kennedy-Lieb-Tasaki (AKLT) Hamiltonian 1 179, 180], can be described as a matrix-product state 
with m = 2. 

For a matrix product state with a particular structure, the elements of the matrices play the role of variational 
parameters. The value of the parameters can, in general, be variationally optimized using a variety of methods. For 
example, the power method (see Sec. 12. 31 . in which the Hamiltonian is repeatedly applied to the variational state, 
ff"|VMp)> imaginary time evolution, e &H \ y/jyip), or iterative diagonalization methods such as the Lanczos method, all 
can be used. The DMRG algorithm is somewhat more complicated to formulate in detail in terms of matrix product 
states; see Ref. 117811 . Expectation values of a set of arbitrary operators residing on different sites can be evaluated 
directly in terms of matrix products if the relevant operators are expressed as m 2 -dimensional operators and are inserted 
into the matrix product at the appropriate points 111 8111 . 

This leads to the question of whether other algorithms based on matrix product states are better than the DMRG 
for particular systems. Recently, a number of such algorithms have been formulated, especially for systems and 
situations in which the DMRG either does not perform well or cannot be applied. In particular, the DMRG does 
not produce translationally invariant states, leading to poor convergence for one-dimensional systems with periodic 
boundary conditions. Verstraete, Porras, and Cirac 1 1 7811 have formulated a variational algorithm that is closely related 
to the DMRG based on a translationally invariant matrix product state. In contrast to the matrix product state generated 
by the DMRG, their state has same degree of entanglement between all sites, including those spanning the boundary. 
The diagonalization of the superblock Hamiltonian used in the DMRG algorithm is replaced by the minimization of 
a more general functional and the sweeping back and forth in the finite system algorithm is replaced by sweeping in 
the positive or negative direction around the ring-like lattice. This algorithm apparently yields a convergence with the 
number of states kept m for the Heisenberg chain with periodic boundary conditions that is comparable to that of the 
DMRG for a chain of the same length with open boundary conditions. However, the computational cost of the new 
algorithm is larger than the DMRG, order m 5 rather than order m 3 because sparseness of the matrix representations is 
lost, and implementations of comparable efficiency to current DMRG codes do not yet exist, so that the utility of the 
method for realistic systems remains to be determined. 



Generalizations of matr ix produc t states have been formulated which incorporate ancilla subspaces so that mixed 
states can be represented Il82lll83il . Such states can be used to represent dissipative quantum or quantum systems 
at finite temperature. Starting with a completely mixed state corresponding to infinite temperature, imaginary time 
evolution with a Trotter-decomposed operator can be carried out to obtain a finite-temperature state 1 182, 183]. Such 
a state can then be further evolved in real time using the methods discussed in Sec. 16.51 yielding the time-dependent 
behavior of finite-temperature or dissipative systems. Since these methods have been proposed quite recently, they still 
await application to realistic systems. 

Ge neralizations of matrix product states that are better suited to two-dimensional systems have also been proposed 
Il04ll . For a two-dimensional square lattice, such a state must represent the entanglement of a site to its four nearest 
neighbors. This can be done by inserting an ancilla space for each bond direction, leading to to a state that is represented 
by a product of rank-4 tensors. (Such tensor product states are related to corner transfer tensor states that have been 
considered for three-dimensional classical systems I127II .') In Ref. 1 1 0411 . such a state is variationally optimized for the 
two-dimensional Heisenberg model with and without frustration using imaginary-time evolution, yielding promising 
variational results for bond dimensions (analogous to the number of states kept, m) of up to 4. However, since the 
algorithm seems to scale with a quite high power of m, its utility for realistic system remains to be explored. 

As we have seen, recent developments in understanding the correspondence between the DMRG and matrix product 
states and in formulating algorithms based on matrix product states and their generalizations have opened up many 
exciting new possibilities in developing new numerical methods for interacting quantum systems. 



6.7. Quantum Information and the DMRG 

As discussed in Sec. 14. 21 the basic approximation used to truncate the Hilbert space of the superblock in the DMRG 
is based on dividing it in two parts given that the superblock is in a pure state (in the simplest case). The accuracy of 
the approximation is intimately connected to the degree of entanglement between the two parts of the bipartite system, 
i.e., between the system block and the environment block. Such entanglement is in turn intimately connected with 
the mutual quantum information in the two subsystems. Recently, progress has been made in utilizing this connection 
between the DMRG and quantum information theory to improve the understanding of the behavior of errors within 
the DMRG, to minimize the error in the DMRG algorithm for a given computational cost, to improve extrapolation 
methods for calculated quantities, and to formulate improved variants of the DMRG algorithm. 

A measure of mutual entanglement is the von Neumann or mutual quantum information entropy S, defined in Eq. 
(1751 . Since this quantity can be written as — w a log w a , its behavior is a measure of the behavior of the eigenvalues 
of the density matrix w a - The von Neumann entropy and related quantities from quantum information theory can 
be used to analyze the behavior of the DMRG algorithms in a quantitative way and to reformulate and optimize the 
algorithms. In particular, Legeza et al. 1 108] have formulated a scheme to vary the number of states kept m dynamically 
according to a criterion based on keeping the truncation error, as measured by the weight of the discarded density- 
matrix eigenvalues, constant. The relationship between the truncation error, the von Neumann entropy, and the size of 
the complete Hilbert space was explored in Refs. 1 109] and 1 184] for the Hubbard model in momentum space and for 
quantum chemistry. The Kholevo bound for the accessible information in a quantum channel, 

X = S{p)-PtgpS(ptw)-(l-pvp)S(p m3P ) , (131) 

where p typ is the density matrix and p tyv the weight of the typical subspace, corresponding to the portion of states 
kept, and p atyp the density matrix of the atypical subspace, corresponding to discarded states, was argued to be a 
good measure of the information loss due to truncation, and was used as a criterion for adjusting m and as a means 
of extrapolating energies 1 1 8-411 . Considerations based on the von Neumann block entropy and the Kullback-Leibler 
entropy |185, 186] have been used to investigate criteria for the optimal ordering of sites and to formulate a mo re 
efficient, dynamical procedure to build up the lattice in the quantum chemistry and momentum-space DMRG 1 109]. 

Another interesting issue is how the behavior of the von Neumann entropy depends on whether or not a system 
is critical. This has been investigated for one-dimensional XY, XXZ, and Heisenberg spin chains in Ref. 1 187]. The 
results are that when the system is not at a quantum critical point, the von Neumann entropy Sl for a particular 
chain length L saturates when L is of the size of the correlation length. For critical systems, Sl is shown to grow 
logarithmically with L, an unbounded behavior, but with a growth that is much slower than the upper bound on the 
entropy, Sl ~ L. The authors argue that the prefactor of this logarithmic divergence depends only on the universality 
class of the critical point. 



The connection between the DMRG and quantum information theory suggests that the DMRG algorithm could 
possibly be applicable to problems in quantum information theory such as how decoherence occurs in detail or what 
the bounds on the information transmitted in a noisy channel are, either through the direct simulation of model systems 
or through the use of insights gained from the analysis of the DMRG algorithm and its generalization using concepts 
from quantum information theory. 

7. CONCLUSION AND OUTLOOK 

In this tutorial, we have discussed a number of closely related numerical methods for investigating strongly interacting 
quantum many-body systems based on numerical diagonalization. The aim is to treat quantum lattice models, systems 
composed of sites, usually on a regular lattice, with a finite number of quantum mechanical states. The size of the 
Hilbert space of such a model grows exponentially with the number of sites, making treatment of sufficiently large 
lattices to scale to the thermodynamic limit difficult. The prototypical models considered here, such as the quantum 
Heisenberg model or the Hubbard model, originate as effective models for interacting electrons in a solid, but quantum 
lattice models can also come about in other physical systems such as atomic nuclei, arrays of trapped atoms, or, in 
general, through the discretization of a quantum field theory. 

All methods discussed here are based on numerical diagonalization, with the NRG and the DMRG having an 
additional variational approximation in which the dimension of the Hilbert space treated numerically is systematically 
reduced. "Exact diagonalization" methods treat the entire Hilbert space of the finite lattice Hamiltonian, which can be 
restricted to particular symmetry sectors. "Complete diagonalization" uses standard algorithms whose cost scales as the 
cube of the dimension to obtain all eigenstates. Since at most Hilbert spaces dimension of order 10 3 , corresponding 
to very small lattice size, can be treated, the method is useful primarily for testing and for benchmarks for other 
methods. "Iterative diagonalization" can calculate extremal eigenstates for much larger Hilbert spaces, up to of order 
10 9 , corresponding to system sizes of up to approximately 40 sites for the simpler quantum lattice models. This is 
done by using an iterative procedure to project the Hamiltonian onto a cleverly chosen subspace. While the method 
is formally variational, convergence to almost machine precision is attained for typical quantum lattice models. The 
methods has also been extended to efficiently calculate dynamical correlation functions, finite-temperature properties, 
and time evolution of a quantum system, properties that, in principle, depend on excited states. The ability of iterative 
diagonalization to provide results that are essentially exact for numerical purposes for a wide variety of properties 
with few restrictions on the type of Hamiltonian or lattice makes it an important standard tool. However, the restriction 
to moderate lattice sizes due to the exponential growth of the Hilbert space of quantum lattice models with system 
size limits its usefulness in extrapolating to the thermodynamic limit. The method is more useful for low-dimensional 
systems simply because the linear dimension is often the relevant parameter in such an extrapolation. 

The numerical renormalization group is an iterative procedure in which a ("complete") numerical diagonalization 
is combined with a truncation of the Hilbert space to lowest energy eigenstates and the addition of degrees of freedom 
by adding lattice sites. The numerical effort is linear in the system size when the number of states kept is held 
fixed. However, the approximation is only well-controlled for one-dimensional quantum lattice models with coupling 
between sites which fall off exponentially with position. Thus, the NRG is limited to systems that can be mapped onto 
such a model; these are primarily single-impurity models such as the Kondo model and the single-impurity Anderson 
model, including its generalization for the dynamical mean-field theory. The NRG is an extremely powerful method 
for such systems and can be used to calculate renormalization group flows, finite-temperature properties and dynamic 
quantities such as the impurity density of states. 

The density matrix renormalization group generalizes the NRG by circumventing problems with the treatment of the 
wavefunction at the boundaries present in the NRG. This is done by embedding a NRG procedure in a larger lattice in 
which the diagonalization is carried out. The reduced density matrix of the renormalized subsystem corresponding to 
one or more states of the system is used to carry out the truncation of the basis in a way that can be shown to be optimal. 
The ability to choose the state or states of the system to which the truncated basis is adapted allows specialization of the 
algorithm for specific purposes, such as calculation of the ground state only or the ground state and low-lying excited 
states or other additional states. The freedom of choosing how to divide the system at each step leads to two classes 
of algorithms: the infinite-system algorithm in which the system grows by two sites at each step, and the finite-system 
algorithm in which the size of the system is kept fixed. While the infinite-system algorithm is useful in building up 
the system at the beginning of the procedure or in scaling directly to the infinite-system limit, use of the finite-system 
procedure leads to more accurate results in almost all cases. 



Crucial to the applicability of the DMRG is how quickly the eigenvalues of the density matrix fall off; this 
determines the accuracy of the variational approximation and is intimately related to the entanglement of the parts 
of the system. For determination of ground-state properties, one-dimensional lattices with open boundary conditions 
can be most accurately treated. It is also favorable for the system to have a gap in one or more symmetry sectors 
of the spectrum, leading to a finite correlation length in the associated correlation function. Critical one-dimensional 
systems have behavior that is less favorable, in particular as the lattice size is increased, and two-dimensional systems 
require exponentially many states to be kept with lattice size in the worst case. Therefore, the DMRG has become an 
extremely powerful standard tool for the calculation of low-energy properties of short-range quantum lattice models 
in one-dimension and on coupled chains. 

Adaptations of the finite-system algorithm to systems with long-range interactions such as the Hubbard model in 
momentum space and quantum chemical systems have been formulated. The momentum-space algorithm becomes 
more accurate as the coupling becomes weaker and its convergence only weakly dependent on dimension; however, 
it is not competitive with the real-space method for one-dimensional and coupled chain systems at moderate to large 
interaction strength. Application to models with longer-range interactions in real space corresponding to shorter-range 
interactions in momentum space should be promising to explore. The quantum chemical DMRG is competitive with 
other quantum chemical methods for small molecules, but more development to improve the performance of the 
algorithm must be done for it to be useful for a wider variety of quantum chemical problems. 

Calculation of dynamical correlation functions using the DMRG can be carried out with adaptations of both methods 
used in iterative diagonalization, the continued fraction (Lanczos basis) method and the correction vector method. The 
latter method provides substantially more accurate information at a given frequency at the price of having to carry out 
a separate calculation at each frequency to be obtained. While some optimization of algorithm is probably possible, 
these techniques seem to be relatively mature, although they remain to be applied widely to one-dimensional quantum 
lattice models. 

Transfer matrices rather than Hamiltonians can also be treated using DMRG methods. For classical systems, 
determining the largest eigenvalue of the transfer matrix allows thermodynamic properties to be calculated for two- 
dimensional discrete systems with an infinite extent in one direction (row transfer matrices) or in both (corner transfer 
matrices) such as Ising or Potts models. Extending the techniques to three dimensions remains a major challenge. 
Treatment of the quantum transfer matrix in the spatial direction for one-dimensional quantum lattice models allows 
thermodynamic quantities and local dynamical quantities to be calculated for infinite lattice sizes. This method is 
accurate in the high-temperature limit and becomes computationally more demanding as the temperature is lowered, 
but is nevertheless applicable to a wide variety of models over a wide temperature range. 

Recent progress in understanding the connection between the DMRG and quantum information theory has opened 
up many new possibilities for the improvement and development of numerical methods. Quantum information can be 
used to optimize the DMRG by providing a criterion for varying the number of states kept to attain a specified error goal 
and to provide insight into optimal ordering of sites for nonlocal models. Efficient calculation of the time evolution of 
quantum states now seems to be possible, a development which has been spurred by insights from quantum information 
theory. Matrix and tensor product states form the basis for possible new algorithms which generalize the DMRG and 
which could potentially attain much better performance for translationally invariant systems and for two-dimensional 
systems and allow the treatment of dissipative systems and better treatment of systems at finite temperature. In addition, 
there is the possibility to apply the DMRG as a numerical technique to improve the understanding of various aspect of 
quantum information theory. 
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